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ABSTRACT 
A  digital  computer  program  which  will  provide  accurate  and  stable- 
solutions  of  linear  and  nonlinear  ordinary  differential  equations  is 
described,  as  are  subprograms  which  simulate  many  of  the  nonlinearities 
occurring  in  feedback  control  systems.   The  program  employs  fourth 
order  Runge-Kutta  numerical  integration  with  automatic  error  checking 
and  interval  modification.   Specialized  coding  sheets  assist  in  the 
preparation  of  input  data,  and  built-in  print  and  graph  output  routines 
provide  a  permanent  record  of  input  equations,,  parameter  values  and 
output  data  for  each  solution.   The  program  is  flexible,,  yet  it  can 
be  used  by  a  person  with  only  a  rudimentary  knowledge  of  FORTRAN 
programming. 
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1.    Introduction. 

The  analysis  of  scientific  problems  frequently  leads  to  a  set  of 
simultaneous  ordinary  differential  equations  which  describe  the  phenomenon 
being  studied.   Analytic  solution  of  these  equations  is  often  impracti- 
cable and  analog  or  digital  simulation  becomes  necessary.   A  digital  com- 
puter program  which  will  solve  ordinary  differential  equations  has  been 
written  in  FORTRAN  60  for  the  CDC  1604  Digital  Computer.   It  includes 
FUNCTION  subprograms  which  cover  many  of  the  nonlinearities  occurring 
in  feedback  control  systems  [7^  •   The  program  has  been  named  SIDES s 
an  acronym  for  Simultaneous  Differential  Equation  Solver. 

Digital  simulation  should  be  considered  as  a  valuable  supplement 
for  the  analog  computer,  but  not  as  a  replacement.   Selection  between 
analog  and  digital  simulation  methods  should  be  based  on  an  evaluation 
of  the  problem  to  be  solved:   the  complexity  of  the  equations,,  the  magni- 
tude of  the  included  parameters  and  the  required  accuracies  of  the  solu- 
tions.  Factors  which  favor  the  selection  of  digital  simulationa  and 
more  particularly,  SIDES,  are: 

a.  The  necessity  for  number  scaling  is  eliminated. 

b.  Solution  accuracies  are  achieved  which  are  unattainable 
with  an  analog  computer. 

c.  Printed  and  graphical  outputs  provide  a  permanent  written 
record  of  input  equations,  initial  parameter  values  and  output  data  for 
each  solution. 

d.  Nonlinear  operations  such  as  squarings  forming  trigono- 
metric and  exponential  functions,  and  multiplication  are  easily  performed. 


e.   Special  functions  such  as  hysteresis  and  time  delay  are 
available  or  can  be  readily  generated  as  ordinary  FORTRAN  SUBROUTINE 
or  FUNCTION  subprograms. 

Some  of  the  limitations  are; 

a.  Inherent  errors  are  due  to  the  sequential  operation  of  the 
digital  computer  in  contrast  to  the  simultaneous  processing  of  signals 
performed  by  the  analog  computer. 

b.  The  effect  of  parameter  variations  on  problem  results  cannot 
be  immediately  observed.   (This  limitation  could  be  partially  eliminated 
by  the  incorporation  of  a  data  display  unit  with  the  CDC  1604). 

Existing  digital  simulation  programs  vary  considerably  with  regard 
to  the  input  language  and  the  integration  techniques  used.   Stein>,  Rose 
and  Parker [5]  have  devised  a  system  which  accepts  as  input  the  descrip- 
tion of  a  "patching"  diagram  prepared  for  the  Electronic  Associates  PACE 
Analog  Computer.   FORTRAN  is  used  as  an  intermediate  in  the  compilation. 
The  DYSAC  program  written  by  Hurley  [2]  is  in  symbolic  machine  language 
(CODAP) .   It  contains  a  supply  of  digitally  simulated  analog  computer 
components  which  are  interconnected  by  the  "patching"  information  sup- 
plied by  the  user.   Integration  is  performed  by  the  fourth  order  Runge- 
Kutta  method  with  a  fixed  integration  interval.   Nordsieck  [4]  has 
developed  a  method  for  automatic  digital  computer  integration  which  is 
equivalent  to  a  reformulation  of  the  method  of  Adams. 

The  SIDES  program  to  be  described  here8  was  designed  with  particular 
emphasis  on  user  convenience  and  solution  accuracy.   Input  and  output 
formats  are  built-in  and  integration  is  performed  by  the  fourth  order 
Runge-Kutta  method  with  automatic  error  checking  and  interval  modification 
as  developed  by  Anderson  [l]  .   A  listing  of  the  SUBROUTINE  and  the  FUNCTION 
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subprograms  which  have  been  developed  for  use  in  conjunction  with  SIDES 
is  given  in  Appendix  B. 


2.    The  SIDES  Input  Routine. 

Calling  program.   SUBROUTINE  SIDES  will  normally  be  called  from  the 
FORTRAN  Library  tape  by  a  calling  program  which  must  be  of  a  prescribed 
form.   To  simplify  the  writing  of  this  program,  a  SIDES  Coding  Sheet  has 
been  prepared  (Appendix  A)  which  reduces  most  of  the  job  of  programming 
to  the  simple  task  of  filling  in  blanks  and  lining  out  undesired  options. 
The  calling  program  deck  must  be  of  the  following  form; 


(Job  ident.) 
(Arbitrary  name) 
(Standard  card) 


..JOB*  JONES  J     BOX  113 
PROGRAM  PROBLEM 

DIMENSION  X(50),  XDOT(30),  C(50) 

C(50)  =1.0  h.     ii 

1  CALL  SIDES  (T,X,XD0T,C)  "     " 

(Equations  written 
in  FORTRAN  language) 
GO  TO  1  (Standard  card) 

END  "     " 

END  " 

This  deck  is  followed  by  a  blank  card  and  the  data  and  option  cards. 

A  typical  program,  together  with  its  output ,  is  included  in 
Appendix  A. 

Equation  statements.   While  equations  for  different  problems  varys 
they  can  be  manipulated  into  the  standard  first  order  form; 


dx 
dt 


i  —  t , ^X- , x_ , • . . , x  , t)     i  —  1  s  /  8 . . .  >,  n 


in  which  the  independent  variable  must  be  designated  by  t  and  the 


dependent  variables  by  x„.   The  number  of  such  equations  is  designated 
by  n  with  n  <^  30.   Since  only  the  x  and  t  variables  can  be  output,  addi- 
tional x.,  n  <  i^.  48,  have  been  made  available.   If  some  quantity  other 
than  one  of  the  variables  is  desired  as  output,  it  must  be  equated  to 
one  of  these  x. . 

A  set  of  constants,  c  ,  1^  i^  49,  has  also  been  made  available  so 
that  parameters  may  be  introduced  or  varied  by  data  cards. 

Any  legal  FORTRAN  statement  or  technique  may  be  used  in  writing 
these  equations,  provided  that  a  transfer  of  control  does  not  prevent  the 
processing  of  all  relevant  statements  during  each  step  of  numerical 
integration.   All  of  the  FORTRAN  Library  functions  are  available  for 
appropriate  inclusion  in  the  equation  statements.   Some  additional 
FUNCTION  subprograms  which  have  been  developed  for  use  by  SIDES  are 
described  below.   These  and  other  subprograms  may  be  appended  to  the 
calling  program  in  the  conventional  manner.   They  may  also  be  added  to 
the  library  tape  if  appropriate. 

Data  and  option  cards.   The  seven  sections  of  input  data  required 
to  specify  a  problem  for  SIDES  are  enumerated  in  Table  I.   All  sections 
are  straight -forward  and  their  purpose  fairly  obvious. 

If  several  solutions  of  the  same  equation  are  required,  it  is 
necessary  to  specify  the  number  of  runs  and  to  supply  a  data  deck  for 
each  run.   For  each  run  after  the  first,  only  the  cards  of  sections  2 
through  7  (Table  I)  are  required.   Information  can  be  retained  between 
runs  by  choosing  the  appropriate  HELD  option.   If  it  is  desired  to  sub- 
divide a  long  solution  into  two  or  more  successive  runs,  the  final  condi- 
tions of  one  run  may  be  used  as  initial  conditions  for  the  next  by  making 

use  of  the  FCS  HELD  command. 
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TABLE  I 


DATA  SECTION 


1.   Problem  description 


2.   Coefficients 


DATA  REQUIRED  FOR  A  PROBLEM 

CONTROL  OPTIONS 
1 


30  first  order  differential 

equations 
1  — *■  9  runs 


Values  zeroed 
Values  set  in 
Values  held 


3.   Initial  conditions 


4.   Time 


Values  zeroed 
Values  set  in 
Initial  or  final  condition 
values  held 


Initial  and  final  time 

given 
Values  held 


5.   Integration  step 
size 


6.   Data  print-out 


7.   Graph  output 


Step  size  given 

Step  size  computed  and 
automatically  varied 
to  maintain  a  pre- 
scribed precision 

Step  size  held 

1 — *■  8  variables  printed 

1 — >  99  steps  of  integration 

between  print-outs 
Data  held 
No  print -out 

1 — *-5  single-curve  graphs 
1 — *-5  curves  on  a  single 

graph 
1— >  99  steps  of  integration 

between  graph  points 
Data  held 
No  graph  output 


3.    The  SIDES  Output  Routine. 

All  job  identification,  together  with  the  run  number  and  a  record 
of  all  input  data,  is  automatically  output  for  each  run.   A  maximum  of 
eight  quantities  may  be  printed-out  during  integration.   An  automatical- 
ly scaled  graph  output  is  also  available,  providing  a  maximum  of  either 
five  single-curve  graphs  or  five  curves  on  a  single  graph. 

A  maximum  of  10,000  integration  incrementss  300  lines  of  print-outs  or 
900  graph  points  (per  curve)  can  be  generated  during  any  one  run.   The  run 
will  be  terminated  if  any  one  of  these  conditions  is  exceeded. 


4.    Error  Indications. 

The  SIDES  program  will  be  terminated  when  all  runs  called  for  are 
completed  or  if  an  error  is  detected  in  the  calling  program  or  in  the 
data  and  option  cards.   SUBROUTINE  SIDES  contains  38  error  stops,  each 
preceded  by  a  diagnostic  print-out.   These  error  stops  refer  to  diffi- 
culties such  as  incorrectly  punched  cardss  missing  cards  and  improper 
option  selection.   All  error  stops  and  print-outs  inherent  to  the  basic 
FORTRAN  Compiler  also  apply. 


5.    The  Integration  Method. 

The  operation  of  SIDES  is  centered  around  the  fourth  order  Runge- 
Kutta  method  of  numerical  integration.   The  choice  of  integration  algo- 
rithm is  necessarily  a  compromise  and  Runge-Kutta  was  chosen  primarily 
because  of  the  ease  both  in  starting  solutions  and  in  changing  the 
integration  increment  during  solutions.   The  use  of  the  method  developed 
by  Anderson  [1\    to  compute  the  truncation  error  has  eliminated  one  of 
the  main  shortcomings  of  Runge-Kutta  as  used  in  other  programs  [6]  „ 
On  this  basis,  the  integration  interval  can  be  automatically  varied  to 
ensure  stability  and  to  keep  the  error  during  each  step  within  prescribed 
bounds,  without  excessive  lengthening  of  the  computing  time.   The  allow- 
able truncation  error  per  unit  time  for  each  of  the  dependent  variables 
was  chosen  to  be: 

-4 
10   (variable  magnitude) 

i  "  total  time 

-9 

The  variation  of  the  step  size  is  bounded  by  arbitrary  limits  of  10 

-3 
and  10  times  the  total  time. 

In  the  event  that  the  user  wishes  to  prescribe  a  fixed  step  sizes 

the  Runge-Kutta  integration  can  be  processed  in  either  one  or  two 

segments,  each  with  a  different  step  size. 


6.    Accuracy  Considerations. 

Since  the  accuracy  of  the  Runge-Kutta  method  of  numerical  integra- 
tion is  a  function  of  the  size  of  the  integration  increment ,  use  of  a 
digital  computer  allows  solution  accuracies  that  in  most  cases  far  ex- 
ceed existing  engineering  demands.   The  following  two  examples  are 
given  to  indicate  the  effect  of  integration  interval  on  solution  accuracy. 
The  interval  was  held  constant  for  each  run. 

Example  1 

2 
Differential  Equation    y'  =  -2xy      (y0  =  1,  x0  =  0) 


Explicit  So! 

Lution 

y  =  (1+x  ) 

1 

Integration 

Interval 

y  at  x  =  5 

Error  at  x  = 

5 

0.03 

3.8xl0'2 

4xl0"12 

0.05 

3.8xl0"2 

-9 
3x10 

0.10 

3.8x10" 

4xl0"8 

2 

Example  2 

Differential 

.  Equation 

y-  =   |2 

y      1+X 

(y0=  1,  x0=  0) 

Explicit  Solution 

y  =  (1+x) 

Integration 

Interval 

y  at  x  =  5 

Error  at  x  =  5 

0.005 

2 
9.4x10   (at 

x  =  2.9) 

-6 
1x10   (at  x=2 

0.03 

8.1xl03 

2xl0"3 

0.10 

8.1xl03 

4x10° 

Milne,  W.  E.  Numerical  solution  of  differential  equations, 
John  Wiley  &  Sons.,  1953:  pp.  20-21. 

2Ibid.,  pp.  74-75. 
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Although  Example  2  is  intentionally  unfavorable  for  solution  by 
Runge-Kutta  integration,  it  is  obvious  that  high  accuracies  can  be 
achieved,  although  at  the  expense  of  additional  computing  time. 

One  of  the  main  advantages  of  the  computed  step  size  integration 
method  is  that  solution  accuracy  and  computing  time  are  not  seriously 
impaired  by  problem  discontinuities.   This  is  an  important  factor  to 
the  engineer  who  must  analyze  systems  containing  switching  devices  or 
step-wise  parameter  fluctuations.   Solution  accuracy  in  SIDES  is  deter- 
mined by  the  estimated  magnitude  (power  of  ten)  of  the  dependent  vari- 
ables, X(j),  1  ^  j  ^  n.   An  examination  of  Example  3  indicates  that 
this  estimated  magnitude  can  often  vary  widely  before  solution  ac- 
curacy suffers. 

Example  3 


where  R  is  of  the  form: 


+5 


-5 


11 


VARIABLE  MAGNITUDE 

C  AT 

ERROR  AT 

C  AT 

ERROR  AT 

SETTING 

T  =  2 
3.2 

T  =  2 
4x10" 10 

T  =  4 

T  =  4 

-3 

- 

-2 

3.2 

5xl0-10 

-2.0 

2xlO-10 

-1 

3.2 

-9 

5x10 

-2.0 

-9 

3x10 

0 

3.2 

lxlO'7 

-2.0 

4xl0"8 

1 

3.2 

2xlO~6 

-2.0 

9xlO~7 

2 

3.2 

7xl0_6 

-2.0 

8xl0_6 

3 

3.2 

7xlO"6 

-2.0 

8xl0_6 

4 

3.2 

7xlO"6 

-2.0 

9xl0_6 
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7.    Special  FUNCTION  Subprograms. 

One  of  the  stated  advantages  of  digital  solution  was  the  relative 
ease  with  which  many  nonlinearities  can  be  simulated.   The  develop- 
ment of  subprograms  to  simulate  functions  such  as  hysteresis  and  time 
delay,  however,  presented  special  problems  in  relation  to  their  use  with 
Runge-Kutta  integration  with  variable  step  size.   Since  FUNCTION  simula- 
tion of  these  nonlinearities  requires  past  information  for  its  proper 
operation,  it  was  necessary  to  modify  the  basic  SIDES  program  so  that 
there  could  be  a  free  interplay  of  information  between  it  and  the 
particular  FUNCTION  subprogram.   Utilizing  the  fact  that  the  COMMON 
statement,  as  applied  to  FORTRAN  60,  will  cause  values  to  be  entered 
in  known  cells,  it  was  possible  to  establish  this  needed  correspondence 
between  the  two  separately  compiled  subprograms. 

FUNCTION  subprograms  have  been  written  to  simulate  time  delay  and 
two  and  three  position  relays  with  hysteresis.   These  subprograms  are 
available  for  use  in  conjunction  with  SUBROUTINE  SIDES.   A  brief  descrip- 
tion of  each  FUNCTION  is  included  in  Appendix  A  and  a  listing  of  each 
FUNCTION  subprogram  appears  in  Appendix  B. 
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8.    Recommendations  and  Conclusions. 

SUBROUTINE  SIDES  has  proven  to  be  an  extremely  easy  program  to  use. 
The  SIDES  Coding  Sheet  enables  the  user  to  prepare  his  program  with  a 
minimum  of  reference  to  additional  explanatory  material.   Due  to  the 
complexity  of  the  input  program,  it  is  believed  that  the  coding  sheet 
is  almost  a  necessity.   Continued  use  of  this  sheet  will  probably  point 
out  needed  changes  to  increase  its  effectiveness.   It  is  already  antici- 
pated that  a  reduction  in  the  data  input  for  successive  runs  is  both 
feasible  and  desirable. 

Although  fourth  order  Runge-Kutta  numerical  integration  with  auto- 
matic error  checking  and  interval  modification  is  accurate,  stable  and 
relatively  rapid,  it  is  recommended  that  other  integration  methods  in- 
cluding predictor-corrector  types,  and  third  and  fifth  order  Runge- 
Kutta  be  evaluated. 

The  FUNCTION  subprograms  developed  for  use  with  SIDES  consist  of  a 
few  of  the  nonlinearities  encountered  in  the  analysis  of  feedback  con- 
trol systems.   Consideration  has  been,  and  should  be  further  given  to 
the  development  of  differentiators,  noise  generators  and  function  genera- 
tors. 

SUBROUTINE  SIDES  is  intended  to  provide  a  simple  but  moderately 
flexible  means  for  the  solution  of  differential  equations.  The  inte- 
gration technique  provides  a  high  degree  of  accuracy  and  stability  and 
the  prepared  coding  forms  simplify  the  programmer's  task.   It  is  be- 
lieved, therefor,  that  SUBROUTINE  SIDES  is  a  valuable  and  essential 
supplement  to  the  analog  computer. 
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SUBROUTINE  SIDES  PROGRAMMING  INSTRUCTIONS 

SUBROUTINE  SIDES  is  intended  to  facilitate  the  solution  of  a  class 
of  problems  involving  linear  and  nonlinear  differential  equations.   It 
was  written  in  FORTRAN  60  for  the  CDC  1604  Digital  Computer.   The  sub- 
routine will  normally  be  called  from  the  FORTRAN  Library  tape  by  a 
program  which  is  being  compiled  and  executed  under  monitor  control. 
This  calling  program  and  the  associated  seven  data  sections  required 
by  SIDES  for  the  solution  of  the  problem  are  supplied  on  80  column 
punched  cards.   Utilization  of  the  SIDES  Coding  Sheets  will  greatly 
simplify  the  programmer's  task  and  serves  as  a  check  that  problem  data 
and  control  commands  are  properly  positioned  and  in  the  correct  order. 
A  description  of  the  functional  sections  of  this  form  and  instructions 
for  their  use  are  given  below.   In  what  follows,  reference  should  be 
made  to  the  copy  of  the  coding  sheets  at  the  end  of  this  appendix. 
These  have  been  filled  out  for  a  typical  problem  and  the  resulting 
computer  output  is  also  appended. 

1.   The  Calling  Program. 

This  section  contains  the  problem  equation  statements  which  will 
probably  cause  the  inexperienced  programmer  the  most  difficulty. 

The  first  card  contains  the  job  identity,  normally  just  the  pro- 
grammer's name  and  his  computer  facility  box  number. 

The  second  card  contains  the  program  name.   Any  combination  of 
seven  letters  or  less  is  acceptable. 

The  next  three  cards  are  standard  and  must  be  punched  as  shown  on 
the  coding  sheet. 
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Following  these  cards  is  the  statement  of  the  differential  equa- 
tions to  be  solved.  These  equations  must  be  written  in  the  form  of  n 
simultaneous  first  order  equations,  where  n  must  not  exceed  30.   Thus, 


XD0T(1)  =  f1(X(l),X(2)f...,X(n),T) 


XDOT(n)  =  fn(X(l),X(2) X(n),T) 

where, 

T       =  the  independent  variable 
X(j)    =  the  jth  dependent  variable 
XDOT(j)  =  the  first  derivative  of  X(j)  with 
respect  to  T 
X(j)  may  be  identified  as  the  output  of  the  jth  integrator  in  the 
corresponding  analog  computer  setup,  while  XDOT(j)  may  be  identified  as 
the  input  to  that  integrator.   Notice,  however,  that  the  sign  reversals 
inherent  in  the  analog  operations  do  not  occur  here. 

It  is  often  desirable  to  define  other  functions  of  the  above  vari- 
ables, either  so  that  the  functions,  f  „,  can  be  built  up  in  more  than 
one  step,  or  so  that  a  quantity  other  than  one  of  the  variables  can  be 
printed  or  graphed.   Although  any  acceptable  FORTRAN  variable  name  may 
be  used  to  define  these  functions,  they  must  be  equated  to  an  X(j), 
n<  j  4  48,  if  they  are  desired  as  output. 

Finally,  if  several  solutions  are  contemplated  with  changes  in  some 
parameters,  the  constants  C(l)  through  C(49)  may  be  introduced.   Their 
values  may  then  be  entered  later  on  data  cards.  The  constant  C(50) 
must  never  be  used  except  as  in  the  calling  program,  above. 
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Block  diagram  manipulations  are  often  useful  in  the  formulation 
of  the  desired  first  order  differential  equations.   The  following 
example  shows  how  a  simple  pole-zero  transfer  function  may  be  reformed 
and  expressed  in  FORTRAN  equation  statements. 


>^C> 


ERROR 


S? 


^Q 


X(2) 


X(l) 


Note 


8+3     m  ,     _3_ 
s+6  s+6 


R  =  C(1)*T 

ERROR  =  R  -   X(2) 

X(2)  =  ERROR  -  X(l)     (permits  print-out  or  graph 

of  C) 
XD0T(1)  =  3.0*(ERR0R  -  2.0*X(1)) 


X(3)  =  ERROR 


(permits  print-out  or  graph 
of  the  ERROR) 


Any  legal  FORTRAN  statement  or  technique  can  be  used  in  the  call- 
ing program,  provided  that  a  transfer  of  control  does  not  prevent  the 
processing  of  all  relevant  statements  during  each  step  of  numerical 
integration.   Improper  transfer  of  control  has  proven  to  be  one  of  the 
most  prevalent  types  of  programming  errors  made  by  users  of  SIDES.   In 
the  variable  step  size  mode  of  integration,  the  equation  statements  are 
frequently  evaluated  with  too  large  a  step  size  and  consequently  too 
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large  a  truncation  error.   These  values  are  then  automatically  reject- 
ed by  the  program  and  the  equations  are  re-evaluated  with  a  smaller 
step  size.   This  process  is  repeated  until  the  truncation  error  is  less 
than  the  prescribed  amount.   Care  must  therefor  be  taken  to  ensure  that 
values  which  have  been  set  after  a  transfer  of  control  on  an  unaccept- 
able pass  are  reset  correctly  on  the  accepted  pass. 

The  COMMON  statement  must  not  be  used  unless  the  programmer  is 
completely  familiar  with  the  SIDES  program.   All  of  the  FORTRAN  Library 
functions  are  available  for  inclusion  in  the  equation  statements.   Addi- 
tional FUNCTION  subprograms  applicable  to  common  control  system  non- 
linearities  have  been  developed  for  use  in  conjunction  with  SIDES 
and  will  be  discussed  later.   If  the  user  desires  to  write  his  own 
FUNCTIONS  or  SUBROUTINES,  they  may  be  added  to  the  calling  program  in 
the  conventional  manner.   Comments  are  permissible  throughout  the  call- 
ing program  if  a  C  is  placed  in  column  one. 

The  calling  program  must  be  terminated  by: 
GO  TO  1 
END 
END 
blank  card 

2.   The  Data  and  Option  Cards. 

C  Values.   One  of  the  options  in  parenthesis  must  be  lined  out. 
The  values  of  the  coefficients,  C(j),  not  simply  zeroed  or  held,  are 
placed,  with  decimal  point,  in  any  order,  on  the  next  three  lines. 
Each  coefficient  is  preceded  by  its  subscript,  j,  right  justified,  in 
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the  appropriate  block.  Unused  lines  are  lined  out  and  the  total 
number  of  lines  (0 — >3)  containing  values  is  noted  on  the  option 
line  for  every  run. 

X  Values.   All  but  one  of  the  options  must  be  lined  out.   On 
any  run  after  the  first,  either  the  initial  or  final  conditions  of 
the  preceding  run  may  be  held  and  automatically  set  as  the  initial 
conditions  of  the  new  run.   The  values  are  entered  in  the  same  manner 
as  the  C  values,  lining  out  the  unused  lines. 

Time.  On  the  initial  run,  the  initial  and  final  problem  times 
are  entered,  with  decimal  point,  and  the  HELD  command  is  lined  out. 
On  subsequent  runs  either  option  may  be  chosen.  The  final  time  has 
an  upper  limit  of  10,000  time  units. 

Integration  Step  Size.  All  but  one  of  the  options  must  be  lined 
out. 

If  the  step  size(s)  are  given,  their  value(s)  are  entered,  with 
decimal  points,  on  the  next  line  with  the  corresponding  termination 
time(s).   The  last  specified  termination  time  must  be  the  same  as  the 
final  time  value  entered  in  the  preceding  section.   Line  out  the  next 
two  lines. 

If  the  step  size  is  to  be  computed,  the  line  for  step  size  data  is 

N 
lined  out  and  the  order  of  magnitude,  N,  (X(j)  ~  10  ,  -9$N^99)  of 

each  dependent  variable,  X(j),  Kj<n,  is  entered,  right  justified, 

following  the  appropriate  variable  subscript  on  the  next  two  lines. 

Unused  variable  subscripts  should  be  lined  out.   If  only  one  dependent 

variable  appears  in  the  problem,  the  second  line  is  lined  out. 
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If  the  HELD  command  is  given,  the  following  three  lines  are  lined 
out. 

Data  Print -out .   All  but  one  of  the  options  must  be  lined  out. 

If  the  GIVEN  command  is  selected,  the  number  of  integration  steps 
between  print-outs  must  be  entered  in  the  appropriate  block.   Although 
a  print-out  at  every  20th  step  will  usually  be  sufficient,  it  should 
be  kept  in  mind  that  the  program  will  automatically  stop  at  an  upper 
limit  of  300  lines  of  print-out.  The  next  line  must  contain  the  sub- 
scripts, right  justified,  of  the  variables  to  be  printed.   For  print- 
out, the  independent  variable,  T,  is  specified  by  the  subscript  50  and 
the  integration  step  size  by  the  subscript  49.   On  the  next  two  lines 
enter  the  data  print-out  column  headings  in  the  appropriate  block  under 
the  corresponding  variable  subscript.   All  three  data  cards  must  be 
present  when  this  option  is  chosen. 

If  the  HELD  or  NONE  commands  are  chosen,  the  following  three  lines 
are  lined  out. 

Graph  Output.   All  but  one  of  the  options  must  be  lined  out. 

If  graph  output  is  needed,  enter  the  desired  number  of  graphs  or 
curves  in  the  appropriate  block.   On  the  next  line  enter  the  number  of 
integration  steps  between  plotting  points.   This  should  be  kept  small 
to  ensure  smoothness,  remembering,  however,  that  the  program  will  auto- 
matically stop  at  an  upper  limit  of  900  plotting  points  per  curve.   On 
the  last  line,  enter  in  the  designated  blocks  the  variable  subscripts 
of  the  X  and  Y  ordinates,  respectively,  for  graphs  (curves)  As  B,  C, 
D,  and  E.   Unused  graph  (curve)  letters  may  be  lined  out.   Since  the 
first  curve  on  the  multi-curve  graph  sets  the  scale,  the  "largest" 
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curve  must  always  be  specified  first. 

If  the  HELD  or  NONE  commands  are  chosen,  the  last  two  lines  are 
lined  out. 

3.    Special  FUNCTION  Subprograms. 

Certain  FUNCTION  subprograms  have  been  developed  for  use  in  con- 
junction with  the  simulation  of  linear  and  nonlinear  feedback  control 
systems  by  SUBROUTINE  SIDES.  These  subprograms  are  called  into  use  by 
the  normal  FORTRAN  techniques  and  may  be  used  with  either  integration 
method  of  SIDES,  The  arguments  may  be  expressed  by  variable  names  or 
numerical  values,  with  decimal  point.  A  brief  description  of  each 
FUNCTION  and  its  arguments  follows. 

FUNCTION  DELAYl  (PELT,  VAR,  T) 

FUNCTION  DELAYl  will  provide  as  output  the  designated  variable's 
value,  delayed  in  time  by  a  prescribed  amount.   The  output  of  the  de- 
lay will  be  the  initial  value  of  the  variable  until  such  time  as  T 
exceeds  DELT.   Since  this  FUNCTION  involves  storage  of  past  informa- 
tion, it  can  only  be  used  once  in  a  given  program.   For  this  reason, 
additional  subprograms,  FUNCTION  DELAY2  and  FUNCTION  DELAY3  are  also 
available.   The  arguments  of  each  of  these  subprograms  are: 

DELT      The  value  of  delay  time.   The  length  of  delay  must 
be  greater  than  the  integration  step  size  and  less 
than  the  sum  of  the  previous  499  time  increments  or 
an  error  print-out  will  be  given  and  the  program 
will  stop. 
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VAR       The  variable  name  of  the  quantity  being  delayed. 
T         The  independent  variable,  time.   It  must  be  represent' 
ed  by  T  in  the  calling  statement. 


FUNCTION  DZONE  (XINPUT,  HALFD,  GAIN) 

FUNCTION  DZONE  simulates  dead  zone  or  free  play. 

OUT 


=  GAIN 


Its  arguments  are: 

XINPUT    The  variable  name  of  the  input  quantity. 

HALFD     The  absolute  value  of  one  half  of  the  dead  zone 

or  band. 
GAIN      The  value  of  the  gain  of  the  nonlinearity. 

FUNCTION  REL2  (XINPUT,  OUTPUT) 

FUNCTION  REL2  simulates  an  ideal  two  position  relay. 

OUT 


OUTPUT 


/ 


IN 
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Its  arguments  are: 

XINPUT    The  variable  name  of  the  input  quantity. 
OUTPUT    The  absolute  value  of  the  output  of  the  relay. 

FUNCTION  REL2H1  (XINPUT.  HALFD,  OUTPUT) 

FUNCTION  REL2H1  simulates  a  two  position  relay  with  hysteresis. 

Since  it  involves  the  storage  of  past  information,  it  can  only  be  used 

once  in  a  given  problem.   For  this  reason,  additional  subprograms, 

FUNCTION  REL2H2  and  FUNCTION  REL2H3  are  also  available. 

OUT 

OUTPUT 


IN 


^-HALFD 


Its  arguments  are: 

XINPUT    The  variable  name  of  the  input  quantity. 

HALFD     The  absolute  value  of  one  half  of  the  hysteresis 

band. 
OUTPUT    The  absolute  value  of  the  output  of  the  relay. 


FUNCTION  REL3  (XINPUT.  OUTPUT) 

FUNCTION  REL3  simulates  an  ideal  three  position  relay. 

OUT 

OUTPUT 


S 


+  IN 
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Its  arguments  are: 

XINPUT    The  variable  name  of  the  input  quantity. 
OUTPUT    The  absolute  value  of  the  output  of  the  relay. 

FUNCTION  REL3D  (XINPUT,  HALFD,  OUTPUT) 

FUNCTION  REL3D  simulates  a  three  position  relay  with  dead  zone, 

OUT 

OUTPUT 


< — r* 


->IN 


HALFD 


Its  arguments  are: 

XINPUT  The  variable  name  of  the  input  quantity. 

HALFD  The  absolute  value  of  one  half  of  the  dead  zone. 

OUTPUT  The  absolute  value  of  the  output  of  the  relay. 


FUNCTION  REL3DH1  (XINPUT,  DROPOUT,  PULLIN,  OUTPUT) 

FUNCTION  REL3DH1  simulates  a  three  position  relay  with  dead  zone 

and  hysteresis.   Since  it  involves  the  storage  of  past  information,  it 

can  only  be  used  once  in  a  given  problem.   For  this  reason,  additional 

subprograms,  FUNCTION  REL3DH2  and  FUNCTION  REL3DH3  are  also  available. 

OUT 


4 

k 

017 

' 

'       i 

\ 

> 

t        * 

k 

\         ^ PULLIN 
^  DROPOUT 

IN 
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Its  arguments  are: 

XINPUT    The  variable  name  of  the  input  quantity. 
DROPOUT   The  absolute  value  of  input  at  which  the  relay 

will  drop  out. 
PULLIN    The  absolute  value  of  input  at  which  the  relay 

will  pull  in. 
OUTPUT    The  absolute  value  of  the  output  of  the  relay. 

FUNCTION  SATAMP  (XINPUT,  CUTOFF,  GAIN) 

FUNCTION  SATAMP  simulates  a  saturating  amplifier. 

OUT 

CUTOFF 


IN 


*i 


J^i_  =  GAIN 


Its  arguments  are: 

XINPUT    The  variable  name  of  the  input  quantity. 
CUTOFF    The  absolute  value  of  the  output  when  it  has 

saturated. 
GAIN      The  value  of  the  linear  gain  of  the  amplifier. 


4.   Sample  Problem. 

The  demonstration  problem  presented  here  illustrates  the  method 

of  data  presentation  required  by  SIDES.   The  completed  SIDES  Coding 

Forms  indicate  the  prescribed  method  of  writing  equation  statements 

and  selecting  program  options.   The  actual  computer  output  appears  next. 

Both  data  print-out  and  single  and  multiple  curve  graphs  are  included  in 

this  output. 
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••JOB*  W.T.  DCOR     BCX  57 

PROGRAM  EXAMPLE 

DIMENSION  X(5G),XD0T<30) ,C(50) 

C(50)=1.0 
1  CALL  SIDES  (T,X,XDOT,C) 

ERROR=C( 1 ) —X ( 1 ) 
C      FCBCK  IS  INCLUDED  tO  DEMONSTRATE  USE  OF  LIBRARY  FUNCTIONS 

FCBCK=ABSF(SINF(X( 1) ))/l 0000.0 

RINPUT=ERROR 

DR0P0UT=C(2)  s 

PULLIN=C(3) 

RELAY=REL3DH2  < RINPUT,C( 2 ) f PULL IN.C (U ) ) 

SIGNAL=RELAY*FDBCK 

XD0T(1)=X(2) 

XD0T(2)=9.0«SIGNAL-U.0«X(2) 

X{3)=ERR0R 

X(M=RELAY 

GO  TO  1 

END 

ENO 


32 


SAMPLE 
OUTPUT 


PROBLEM 
DATA  FCR 


RUN 


W.  T 
NUMBER 


1 


DOOR  BOX  57 


TIME 

OUTPUT 

OUTPUT 

ERROR 

RELAY 

POSITION 

VELOCITY 

OUTPUT 

.OOOOOE+OO 

-.500C0E+00 

•1OOO0E+0O 

•15000E+01 

.15000E+01 

.1U270E+00 

-.37Uii8E  +  00 

•152UUE+01 

.137U5E+01 

, .15000E+C1 

.34270E+00 

.U5759E-01 

•25U35E+01 

.95U2UE+00 

.1500CE+01 

.5U270E+00 

.60629E+00 

•3001UE+01 

.39371E+00 

.15C00E+01 

.59006E+00 

.750COE+00 

.30659E+01 

.25000E+00 

.0OCO0E+0O 

.65054E+00 

.91U72E+00 

..24071E+01 

.85283E-01 

.OOOOOE+CC 

.85054E+00 

•12U61E+01 

.10817E+01 

-.2U611E+00 

.OOOOCE+OO 

•10505E+01 

•13950E+01 

•48615E+00 

-.39503E+00 

.CO0C0E+0O 

•12505E+01 

•1U620E+01 

•21857E+CO 

-.46198E+00 

.OOOOOE+OO 

.1U505E+01 

•14921E+01 

•98331E-01 

-.M9208E+00 

.OOOOOE+CO 

.15U77E+01 

•150C0E+01 

•66736E-01 

-.50000E+00 

.OOOOOE+CO 

.15U90E+01 

•150C1E+01 

•U8739E-01 

-.50008E+00 

-.15000E+C1 

.17282E+G1 

•13333E+01 

-.17029E+01 

-.333.30E+00 

-.1500GE+01 

•17733E+01 

•125COE+01 

-.19791E+01 

-.25000E+00 

,  .COOOOE+OO 

•18238E+01 

•  115<35E+01 

-.161 72E+01 

-.15953E+00 

.OOOOOE+OO 

.20238E+01 

•93691E+00 

-.72653E+00 

.63092E-01 

.OOOOOE+CC 

.22238E+01 

•83690E+00 

-.32636E+00 

.16310E+00 

•COOOOE+OO 

.2U238E+01 

•79198E+00 

-.1U655E+CO 

•20802E+00 

.COCOOE+OO 

.26238E+01 

•77181E+00 

-.65763E-01 

.22819E+00 

.OOOOOE+OO 

•28238E+01 

•76277E+00 

-.29U63E-01 

.23723E+00 

.OOOOOE+CO 

•30238E+01 

•75872E+00 

-.13153E-C1 

.2U128E+00 

.OOOOOE+OO 

.32238E+01 

•75692E+00 

-.5825CE-C2 

.24308E+00 

.OOOOCE+CO 

•34238E+01 

•75613E+00 

-.25323E-02 

.2U387E+00 

.COOOOE+OO 

•36238E+01 

•75579E+00 

-.10529E-C2 

.2U1421E  +  00 

.COOOOE+OO 

.38238E+01 

•75566E+00 

-.3881 1E-03 

.2443UE+00 

.OOOOOE+CC 

.U0238E+01 

•75561E+00 

-.89U23E~0U 

.2U439E+00 

.OOOOCE+OO 

•42238E+01 

•75561E+00 

.I4U7  82E-OU 

.24U39E+00 

.COOOOE+OO 

•UU238E+01 

•75563E+00 

•10509E-03 

.24437E+00 

.OOOOOE+OO 

.U6238E+01 

•75565E+00 

•13218E-03 

•2U435E+00 

•OOOOOE+CC 

.U8238E+01 

•755d8E+00 

•1UU36E-03 

•24432E+00 

•OOOOOE+OO 

NORMAL  STOP  AT  FINAL  TIME 


GRAPH  TITLEO 


SAMPLE  PROBLEM 


W.T.  DOOR*  BOX  57   RUN  1   GRAPH  A 


GRAPH  TITLED  ..   SAMPLE  PROBLEM 


W.T.  DOOR  BOX  57   RUN  1   GRAPH  E 


Printout  for  second  run  not  Included. 
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SAMPLE    PROBLEM 

2    RUNS    ARE    CALLEC    FOR. 


W. T.  DOOR  BOX  57 


INPUT  OATA  FOR  RUN  NUMBER  1 


THE  NON-ZERO  DATA  COEFFICIENTS  ARE  AS  FOLLOWSt 


C( 

1) 

s 

.1COE+01 

C( 

2) 

3 

•250E+00 

C( 

3) 

■ 

.500E-K30 

C( 

U) 

S 

•150E+01 

THE  NON-ZERO  INITIAL  CONDITIONS  ARE  AS  FOLLOWS* 


XO( 
XO( 


1)  = 

2)  * 


-.500E+00 
•lOOE+OO 


THE  TIMING  DATA  ARE  AS  FOLLOWS, 


INITIAL  TIME 
FINAL  TIME 


.OOOOOE+OO 
•5OC00E+O1 


THE  STEP  SIZE 
THE  ESTIMATED 
THE  VARIABLES 


IS  COMPUTED,  BASED  ON 
ORDER  OF  MAGNITUDE  OF 
AS  LISTED  BELOW. 


PRINT 
20 


GRAPH 


X(  1)  * 
X(  2)  * 

l.OE  Q 
l.OE-3 

K 

SUMMARY 

m 

INCREMENTS 
E  VARIABLES 

BETWEEN 
PRINTED 

PRINTOUTS. 
ARE, 

X(50) 
X(  1) 
X(  2) 
XC  3) 
X(  4) 

TIME 

OUTPUT 

OUTPUT 

ERROR 

RELAY 

POSITION 
VELOCITY 

OUTPUT 

SUMMARY 

• 

2  SINGLE  GRAPHS. 

PLOT  EVERY  PCINT. 

GRAPH  A  IS  X(50)  VS.  X(  1) 
GRAPH  B  IS  X(50)  VS.  X(  3) 

•  X(50)  REPRESENTS  TIME. 


- i  l~M|P  TIT^HIIIMP" 


■  *m 


SAMPLE    PROBLEM  W.T.    DOOR    BOX    57 

OUTPUT    DATA    FOR    RUN   NUMBER    1 


TIME 

OUTPUT 

OUTPUT 

ERROR 

.  RELAY 

POSITION 

VELOCITY 

1 

OUTPUT 

.OOOOOE+OO 

-.500C0E+00 

.10000E+0G 

.15000E+01 

.15000E+01 

.1U270E+00 

-.37U148E  +  00 

.152U4E+01 

.13745E+01 

.15000E+C1 

.34270E+00 

.U5759E-01 

.25U35E+01 

.95U2UE+00 

.1500CE+01 

.5U270E+00 

.60629E+00 

.3001UE+01 

.39371E+00 

.15C0GE+01 

.59006E+00 

.750C0E+00 

.30659E+01 

•25000E+00 

.OOOOOE+00 

.6505UE+00 

•91U72E+00 

.2U071E+01 

.85283E-01 

.OOOOOE+OO 

.85054E+00 

.12U61E+01 

.  108 17E+01 

-.2U611E+00 

•OOOOOE+00 

.10505E+01 

•13950E+01 

.M8615E+00 

-.39503E+00 

.COOOOE+OO 

.12505E+01 

•1U620E+01 

.21857E+CC 

-.46198E+00 

.OOOOOE+OO 

.14505E+01 

.14921E+01 

.98331E-01 

-.U9208E+00 

.OOOOOE+CC 

.15U77E+01 

.150C0E+01 

.66736E-01 

-.50000E+00 

.OOOOOE+00 

.15U90E+01 

.150C1E+01 

.48739E-01 

-.50008E+00 

-.15000E+C1 

.17282E+G1 

•13333E+01 

-.17029E+01 

-.33330E+00 

-.1500CE+01 

•17733E+01 

.125C0E+01 

-.19791E+01 

-.25000E+00 

.COCOOE+00 

.18238E+01 

.11595E+01 

-.16172E+01 

-.15953E+00 

.OOOOOE+OO 

.20238E+01 

.93691E+00 

-.72653E-I-00 

.63092E-01 

.OOOOOE+CC 

.22238E+01 

•83690E+00 

-.32636E+00 

.16310E+00 

.OOOOOE+00 

.2U238E+01 

.79198E+00 

-.14655E+C0 

.20802E+00 

.COOOOE+OO 

.26238E+01 

.77181E+00 

-.65763E-01 

.22819E+00 

.OOOOOE+00 

.28238E+01 

.76277E+00 

-.29U63E-01 

.23723E+00 

•OOOOOE+00 

.30238E+01 

.75872E+00 

-.13153E-C1 

.2U128E+00 

.OOOOOE+00 

•32238E+01 

.75692E+00 

-.5825CE-C2 

.2U308E+00 

.OOOOCE+CO 

.34238E+01 

.75613E+00' 

-.25323E-02 

.2U387E+00 

.COOOOE+OO 

.36238E+01 

.75579E+00 

-.10529E-C2 

.2UH21E+00 

.COOOOE+OO 

.38238E+01 

.75566E+00 

-.3881 1E-03 

.24M3UE+00 

•OOOOOE+OC 

•U0238E+01 

.75561E+00 

-.89U23E-0U 

.2U439E+00 

.OOOOCE+00 

•42238E+G1 

•75561E+00 

•UU782E-0U 

.'24U39E+00 

.COOOOE+OO 

•44238E+01 

.75563E+00 

.10509E-03 

•  24J437E  +  00 

.COOOOE+OO 

.U6238E+01 

.75565E+00 

.13218E-03 

.2U435E+00 

.OOOOOE+CC 

.U8238E+01 

•75568E+00 

.1U436E-03 

.24U32E+00 

*00000E+00 

NORMAL    STOP    AT    FINAL    TIME 


GRAPH    TITLED    ..       SAMPLE    PROBLEM 


W.T.    DOOR    BOX    57      RUN    1       GRAPH    A^ 


GRAPH   TITLED    ..       SAMPLE    PROBLEM 


W.T.    DOOR    BOX    57      RUN    1       GRAPH    B 


Printout  for  second  run  not  Included. 
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APPENDIX  B 


38 


LIST!  NG  0  F  S  UB R  GU  T I NE  SI  D  E  S 


SUBROUTINE  SIDES  (T,  X,  XDCf,  C) 

DIMENSION  X(50),  XDOT130),  CC50),  II;TcLi10»,  JTITLE  Mo),  XRI30), 

1  ip(IO),    ir. i20)  •   xo{30),   PRHoi,   G'uioi,   xic500j»   crui, 

2  YH9CC),    X2I90G),     Y2J90C5,     X3I900),     Y3S9CC),     X4i9CC), 

3  Y4 I  900  )  ,     X  5  J  90  0  » ,     Y5  f  900 } ,     A { 3 1) ,     A K *  4 , 30  )  , F I  6 )  , 

4  XS(30)f     XRRI30  5,    XAI30),     XB130),    KPI30),     ICHECM30}, 

5  JP ( 30  )  , V  (1  500  J  t W i) 500  I , MY  5  30), D i 40  ) 

COMMON  XI ,X2,X3,X4,X5,Y1  ,  Y2,  Y3,  Y«* ,  Yb  ,  ICHECK,  I TI  ILE, J  I  I  TLE  , 

1  AK,F  , A , X  S , XR , XRR , X A , X  B , KP , JP , XO,  1 P , P R , V , W *  FY , C 

1000  FORMAT  (1H1) 

1010  FORMAT  (/} 

1020  FORMAT  (// ) 

1030  FORMAT  (///) 

1040  FORMAT  (1GA8) 

1050  FORMAT  {  12 , A6, A5 , II , A8) 

1060  FORMAT  (36H  EQUATIONS/ RUNS  DATA  CARD  INCORRECT.) 

1070  FORMAT  [35H  NUMBER  OF  EQUATIONS  NOT  SPECIFIED.) 

1071  FORMAT  (41H  MAXIMUM  NUMBER  OF  EQUATIONS  Al LHWED  IS  30.) 
1080  FORMAT  (30H  NUMBER  OF  RUNS  NOT  SPECIFIED.) 

1090  FORMAT  I5X.22HCNE  RUN  IS  CALLED  FOR./, 

1100  FORMAT  (5X,I1,21H  RUNS  ARE  CALLED  FOR./, 

1110  FORMAT  I31H  MORE  THAN  NINE  RUNS  REQUESTED.) 

1120  FORMAT  ( AR, A3, A6,2A4,4A8 , A7,  11 , A8) 

1121  FORMAT  (38H  ON  C  VALUE  CARD  CROSS  OUT  ONE  OPTION.) 
1130  FORMAT  (35H  C  VALUE  CARD  INCORRECT  OR  MISSING.) 
1140  FORMAT  (8{ I3.F7.2) ) 

1150  FORMAT  (27H  INPUT  DATA  FOR  RUN  NUMBER  ,I1,//',47H  THE  NON-ZERO  HAT 

1A  COEFFICIENTS  ARE  AS  FOLLOWS,/) 

1160  FORMAT  I10X, 2HC ( , I  2, 4H )  =  ,£8.3) 

1170  FORMAT  (45H  MAXIMUM  ALLOWED  COEFFICIENT  SUBSCRIPT  IS  49.) 

1180  FORMAT  (5X,4HN0NE) 

1190  FORMAT  i A8,A3, A6, A4, A7,A3, A7,3A8, A2, II • A8) 

1191  FORMAT  {4UH  ON  THE  X  VALUE  CARD  SELECT  ONE  OPTION  ONLY.) 
1200  FORMAT  (35H  X  VALUE  CARD  MISSING  OR  INCORRECT.) 

1210  FORMAT  (48H  MAXIMUM  ALLOWED  STATE  VARIABLE  SUBSCRIPT  IS  30.) 

1220  FORMAT  (48H  THE  NON-ZERO  INITIAL  CONDITIONS  ARE  AS  FOLLOWS,/) 

1230  FORMAT  II  OX, 3HXC I , 12 ,4H)  =  ,118.3) 

1240  FORMAT  ( A8, A3, A8 , A4, F7.3 , A8, A5, F7. 3, A8, A3, A4 , A8) 

1250  FORMAT  (50H  ON  THE  TIME  DATA  CARL  ONE  OPTION  NOT  CROSSED  CO f . I 

1260  FORMAT  (29H  ERROR  ON  THE  TIME  DATA  CARD.) 

1261  FORMAT  ( A8, A3, A5  , A8, A7,3 A8, A6, A4  ,A8) 

1270  FORMAT  (51H  SELECT  ONE  OPTION  ONLY  ON  THE  STEP  SIZE  DATA  CARL.) 

1280  FORMAT  ( I  1 , I  2 , 12  , 1 4 (  13, 1  2  )  ) 

1290  FORMAT  (77H  STEP  SIZE  MODE  AND  GIVEN  DATA  ARE  NOT  CONSISTENT  OR  TH 

1ERE  IS  A  FORNAT  ERROR.) 

1300  FORMAT  (48H  STEP  SIZE  OPTION  NOT  INDICATED  OR  IS  INCORRECT.) 

1310  FORMAT  (2  I  I  1  ,2A6,F7.3, A6, A7, F7. 3 ) ) 

1320  FORMAT  (51H  FOR  GIVEN  STEP  SIZE  OPTION,  VALUE  OF  CT  NOT  GIVEN. 9 

1330  FORMAT  (54H  FOR  GIVEN  STEP  SIZE  OPTION,  VALUE  OF  IF  MOST  BE  GIVEN) 

1340  FORMAT  (32H  THE  TIMING  DATA  ARE  AS  FOLLOWS,,// 

1  5X,15HINITIAL  TIME  =  ,E11.5,/ 

2  5X,15HFINAL  TIME    =  ,E11.5) 

1350  FORMAT  {5X.35HTFE  STEP  SIZE  IS  COMPUTED,  BASED  ON,/ 

1  5X,36FTHt  ESTIMATED  ORDER  OF  MAGNITUDE  OF  / 

2  5X,30FTHE  VARIABLES  AS  LISTED  BELOW.//) 
1360  FORMAT  MOX  ,2HX  (  ,  I  2,  4H  )  =  ,4Hl.0E,I2) 

1370  FORMAT  (5X.15HSTEP  SIZE     ■  , E 1 1 . 5 , 1 0H  FROM  T  =  ,E11.5,8N  TC  T  = 

1  ,E11.5) 

1380  FORMAT  ( A8, A3,3A8, A7, 12, 2A8, Al , A4, A6, A4, A4) 

1390  FORMAT  (39H  PRINTOUT  CARD  IS  MISSING  OR  INCORRECT.) 
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1U00  FORMAT  C42H  CN  PRINTOUT  CARD,  SELECT  ONE  OPTION  ONLV„S 

1410  FORMAT  (16,7110) 

1420  FORMAT  (II  ,8(A8,I2)) 

1430  FORMAT  (A8,A3,I 1 ,2A8,A3,  I1,3A8,A5,A4,A6»A4,A2> 

1440  FORMAT  (38H  CRAPH  DATA  CARD  MISSING  OR  INCORRECT.) 

1450  FORMAT  (38H  CN  GRAPH  CARD  CHOOSE  ONE  OPTION  ONLY.) 

1460  FORMAT  ( A8 , A3 , 12 , A8 ) 

1470  FORMAT  (8 ( A2  ,  12 , 12 , 12, 12 ) ) 

1480  FORMAT  {44H  MAXIMUM  NUMBER  OF  CURVES  OR  GRAPHS  IS  FIVE.) 

1490  FORMAT  (48H  SPACING  BETWEEN  GRAPH  POINTS  MUST  BE  SPECIFIED.) 

1500  FORMAT  (46H  NUMBER  OF  CURVES  OR  GRAPHS  MUST  BE  SPECIFIED.) 

1510  FORMAT  (52H  PLOT  EVERY  --TH  POINT  CARO  IS  MISSING  OR  INCORRECT.) 

1520  FORMAT  (17H  PRINT  SUMMARY ,/) 

1530  FORMAT  (5X.11HN0  PRINTOUT) 

1540  FORMAT  (5X,I2,3CH  INCREMENTS  BETWEEN  PRINTOUTS./ 
1         5X,    26HTHE  VARIABLES  PRINTED  ARE,/) 

1541  FORMAT  (5X.18HPRINT  EVERY  POINT./ 

1         5X,    26HTHE  VARIABLES  PRINTED  ARE,/) 

1542  FORMAT  (5X,I2,3CH  INCREMENTS  BETWEEN  PRINTOUTS./ 
1         5X,    2UHTHE  VARIABLE  PRINTED  IS,/) 

1543  FORMAT  (5X,18HPRINT  EVERY  POINT./ 

1         5X,    24HTHE  VARIABLE  PRINTED  IS,/) 
1550  FORMAT   ( 1  OX ,2HX { , 12 . 1 H)  , 5X, A8, 1 X, A8 ) 
1560  FORMAT  J/,5X,24H*  X { 50 )  REPRESENTS  TIME.) 
1570  FORMAT  (17H  GRAPH  SUMMARY—,/) 
1580  FORMAT  (5X,1CHN0  GRAPHS.) 

1590  FORMAT  <6X,I1,1SH  SINGLE  GRAPHS.) 

1591  FORMAT  (6X,17H0NE  SINGLE  GRAPH.) 

1600  FORMAT  (6X.4CH0NE  MULTIPLE  GRAPH  LABELED  GRAPH  M  WITH  ,11,  16H  CUR 
IVES  PLOTTED.) 

1601  FORMAT  (58H0NE  MULTIPLE  GRAPH  LABELED  GRAPH  M  WITH  ONE  CURVE  PLCTT 
1ED.  ) 

1610  FORMAT  (5X,I2,27H  INCREMENTS  BETWEEN  POINTS./) 

1611  FORMAT  (6X,17HPL0T  EVERY  POINT.) 

1620  FORMAT  ( 1 0X, 1 3HGRAPH  A  IS  X(,I2,8H)  VS.  X(,I2,lH)) 
1630  FORMAT  { 1 0X, 1 3HGRAPH  B  IS  X(,I2,8H)  VS.  X(,I2,1H)) 
1640  FORMAT  ( 10X, 1 3HGRAPH  C  IS  X(,I2,8H)  VS.  X(,I2,1H)) 
1650  FORMAT  { 10X, 13HGRAPH  D  IS  X(,I2,8H)  VS.  X(,I2,1H)) 
1660  FORMAT  ( 1  OX, 13HGRAPH  E  IS  X(,I2,8H)  VS.  XC , I  2, 1H) ) 
1670  FORMAT  { 10X, 1 3HCURVE  A  IS  X(,I2,8H)  VS.  X(,I2,lH)) 
1680  FORMAT  { 10X, 1 3HCURVE  B  IS  X(,I2,8H)  VS.  X(,I2,lH)) 
1690  FORMAT  ( 1 0X,  13HCURVE  C  IS  X(,I2,8H)  VS.  X(,I2,1H)) 
1700  FORMAT  { 1  OX, 1 3HCURVE  D  IS  X(,I2,8H)  VS.  Xl,I2,lH)) 
1710  FORMAT  t 10X, 13HCURVE  E  IS  X(,I2,8H)  VS.  X(,I2,lH)J 
1730  FORMAT  (4X,1CA8) 

1740  FORMAT  (5X, 27HOUTPUT  DATA  FOR  RUN  NUMBER  ,11//////) 

1741  FORMAT  {  1X,5HPAGE  ,I1,1X,30H0F  OUTPUT  DATA  FOR  RUN  NUMBER  ,11,/ 
1        ////////) 

1750  FORMAT  < 3X, 7 ( A8, 4X ) , A8 ) 

1760  FORMAT  (10( IX, El  1.5)  ) 

1770  FORMAT  (//24F  STOP  AT  300  PRINT  LINES) 

1780  FORMAT  I//19H  STOP  AT  T  =  10,000) 

1790  FORMAT  (//26H  STOP  AT  10,000  INCREMENTS) 

1800  FORMAT  {//25F  STOP  AT  900  GRAPH  POINTS) 

1810  FORMAT  (//26H  NORMAL  STOP  AT  FINAL  TIME) 

1820  FORMAT  <//llF  STOP  AT  X(,I2,10H)  =  10,000) 

1830  FORMAT  (30H  ALL  RUNS  HAVE  BEEN  COMPLETED.) 

IND=C(50) 

GO  TO  (2000, 361C, 3680, 5020, 3523, 3530) ,IND 

2000  DO  2030  J=1,  10 
2030      ITITLE(J)=8H 

DO  2040  J=l,16 
2040      JTITLE(J)=8H 
C      RK  I  COEFFICIENTS 

CT(1  )=0.0 

CT(2)=0.5 

CT(3)=0.5 

CT(4)  =  1  .0 
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NRC  =  0 

PRINT  1000 
C      READ  TITLE  CARD 

READ  1040,  (  ITITLEU  ),  1=1,5) 
C      READ  NUMBER  CF  RUNS  /  EQUATIONS  CARD 

READ  1050,  \E,  ICHECK ( 1 )  ,  ICHECK { 2) , NR f  ICHECK ( 3 ) 

ICHECK(4)=8H  ECOATIC 

ICHECK(5)=8H  RUNS 

IF  {  ICHECK(4)-ICHECK(1 ))  2060,2050,2060 
2050  IF  UCHECK«5)-ICHECK(3))  2060, 2C70, 2060 
2060  PRINT  1060 

STOP  1 
2070  IF  (NE)  2080,2080,2090 
2080  PRINT  1070 

STOP  2 

2090  IF  (NE-30)  2092,2092,209  1 

2091  PRINT  1071 
STOP  3 

2092  PRINT  1730,  ( IT  ITLE (  I ) , I = 1 , 5 ) 
IF  (NR-1)  2100,2  110,2120 

2100  PRINT  1080 

STOP  4 
2110  PRINT  1090 

GO  TO  2130 
2120  PRINT  1100, NR 
2130  PRINT  1030 

IF  (NR-9)  2150, 2150, 2140 
2140  PRINT  1110 

STOP  5 
C      ENTRY  POINT  FOR  ADDITICiMAL  RUNS 

2150  NRC=NRC*1 
INDI=1 
NPAGE=1 
MZ=1 

C      READ  C  VALUE  CARD 

READ  1120, (I  CHECK ( I ) , I =1  , 10 ) , NC, ICHECK!  11  ) 

I  CHECK! 15)=6KZERCED 

ICHECM  16)=4HHELD 

I  CHECK  {  17)=ICHECK(  15)-UCHECK(  16) 

I  CHECK ( 18)=ICHECK(3)+ICHECK(5) 

IF  ( ICHECKl 1 7 )-I CHECK ( 18 ) )2 1 52, 2 151 , 2 1 52 

2151  PRINT  1121 
STOP  6 

2152  IF  ( ICHECK(3)-ICHECK(15) )  2160,2180,2160 
2160  IF  l ICHECK{5)-ICHECK(16) )  2170,2200,2170 
2170  PRINT  1 130 

STOP  7 
2180  DO  2190  J=1,»49 
2190      CU)=0. 
2200  IF  (NC)  2220,2260,2220 
C      READ  C  VALUES 
2220  DO  2250  J=1,NC 

READ  1140  i JP(K) ,XA{ K),K=1,8) 
DO  2250  K=1,8 

IF  UP(K)-49)  2240,2240,2230 
2230        PRINT  117C 

STOP  8 
2240        ITK=JP(K) 
2250        C(ITK)=XA(K) 
2260  PRINT  1150,  NRC 
K  =  0 
DC  2280  J=1,M9 

IF  (C( J) )  2270,2280,2270 
2270      PRINT  1160, J, C(J) 

K=K  +  1 
2280      CONTINUE 

IF(K)  2300,2290,2300 

PRINT  1180 

PRINT  1020 

READ  X  VALUE  CARD 
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RE AC    1  190, C  I  CHECK (J) ,J  =  1  ,11  ),NI , I  CHECK! 12) 

ICHECKt 13)=6FZERGED 

ICHECKt  14)=7FIC    HELD 

ICHECKt  15)=7hFC    HELD 

I  CHECK ( 16)=ICHECK( 13)  + ICHECKt 14) 

I  CHECK ( 17 )  =  I CHECK ( 13)+ ICHECKt 15) 

I  CHECK <  18)  =  ICHECKt 14 )+ [CHECK { 15) 

ICHECKt  19  )  =  ICHECKt  3)  +  ICHECM5) 

ICHECK«20)=ICHECK(3)+ICHECK{7) 

ICHECKt 21 )=ICHECK(5) + ICHECKt 7) 

IF    <  ICHFCKt  16)-ICHECK{ 19) )     2301,2303,2301 

2301  IF    t  ICHECKt 1 7 )- ICHECKt 20 ) )    2302,2303,2302 

2302  IF    t ICHECKt 1 8 )- ICHECK t 21 ) )    2304 , 2303, 230* 

2303  PRINT    1 191 
STOP    9 

2304  IF    t  ICHECKt3)-ICHECKt13)  )    2330,2310,2330 
2310    00    2320    J=1,NE 

2320      XCtJ)=0. 

GO  TO  2380 
2330  IF  (ICHECKt 1  4  )- ICHECK ( 5)  )  2340,2380,2340 
2340  IF  (  ICHECKt 1 5 )- ICHECKt 7)  )  2350,2360,2350 
2350  PRINT  1200 

STOP  10 
2360  MZ  =  2 

DC  2370  J=1,NE 
2370      XOtJ)=XtJ) 
2380  IF  (NI)  2410,2430,2410 
C      READ  X  VALUES 

2410  DC  2420  J=1,M 

READ  114C  UPtK)  ,XAt  K),K=1,8) 
DO  2420  K=1,8 

ITK=JP{K) 

IFtITK-30)  2420,2420,2411 

2411  PRINT  1210 
STOP  12 

2420        XOC ITK)=XA(K) 
2430  PRINT  1220 

K=0 

DO  2450  J=1,NE 

IF(XOtJ))  2440,2450,  2440 
2440      PRINT  1230,  J,  XOtJ) 

K  =  K+1 
2450      CONTINUE 

IFtK)  2470,2460,2470 
2460  PRINT  1180 
2470  PRINT  1020 
C      REAC  TIME  DATA  CARD 

READ  1240, { ICHECKt J) , J  =  1  ,4), TO, ICHECK t 5 ), ICHECK ( 6 ), TF ,( ICHECK 1 1 ),  I 
1=7,10) 

ICHECKt  11)=8HSTART  TI 

ICHECKt 12)=4FHELD 

ICHECKt 13)=ICHECKt 11 )♦ ICHECKt 12) 

ICHECKt  14)=ICHECM3)+ICHECK{9) 

IF    (  ICHECK113)-ICHECKM4)  )    2490,2480,2490 
2480    PRINT    1250 

STOP  13 

2490  IF  t  ICHECKt 1 2 )- ICHECKt 9)  )  2500,2491,2500 

2491  TF=F(1) 
TC=F(2) 

GO  TO  2530 
2500  IF  t  ICHECKt 11)-ICHECKt3)  )  2520,2521,2520 

2520  PRINT  1260 
STOP  14 

2521  F(1)=TF 
F(2)=T0 

C      REAC  STEP  SIZE  CARD 
2530  READ  1261,  t  ICHECK t I ), 1= 1 , 1 1 ) 
ICHECKt 12)=5HGIVEN 
ICHECKt 13)=8HC0KPUTED 
ICHECKt 14)=4HHELD 
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I  CHECK { 15 )=I CHECK* 12 ) + ICHECK C 1 3 ) 

I  CHECK ( 16)=ICHECK( 1 2 )♦ ICHECK ( 1 4 ) 

I  CHECK  {  17)  =  ICHECK {  1  3  )  +  IC HECK { 14 ) 

I  CHECK (  18)=ICHECK(3)+ICHECK(6) 

ICHECK ( 19)=ICHECK(3»+ICHECK( 10) 

ICHECK(20)=ICHECK(6)+ICHECK(10) 

IF  (  ICHECK(15)-ICHECK(18)  )  2540,  2560  ,  2540 
2540  IF  (  ICHECKC 16)-ICHECK( 19) )  2550,2560,2550 
2550  IF  (  ICHECKI 1  7  )-ICHECK ( 20 ) )  2570,2560,2570 
2560  PRINT  1270 

STOP  15 

IF  (ICHECKi 13)-ICHECK(6) )  2610,2580,2610 

READ  COMPUTEC  STEP  SIZE  CARD 

READ  1280,  NC,  ( NH, KP ( I ) , I =1 , NE, 2 ) 

IF  (ND-2)  2590,2600,2590 
2590  PRINT  1290 

STOP  16 

IF  (NE-1)  2601,2720,2601 

READ  COMPUTEC  STEP  SIZE  CARD 

READ  1280,  NC,  ( NH,KP (I ) , 1=2, NE, 2) 

IF  (ND-2)  2590,2720,2590 
2610  IF  ( ICHECK(1M)-ICHECK( 10) )  2620,2720,2620 
2620  IF  C  ICHECK! 1 2 )- ICHECKI 3)  )  2630,2640,2630 
2630  PRINT  1300 

STOP  17 
;      READ  GIVEN  STEP  SIZE  CARD 
2640  READ  1310 , ND, NH , NH,F { 3 ) , NH, NH ,F ( 4 ) ,NH , NH, NH, F { 5 ) ,NH,NH,F(  6 ) 

X(49)=F(3) 

IF  (ND-1)  2650,2660,2650 
2650  PRINT  1290 

STOP  18 
2660  IF  (F(3))  2680,2670,2680 
2670  PRINT  1320 

STOP  19 
2680  NCT=1 

IF  (F(4)  -TF)  2690,2  720,2690 
2690  IF  (F(5))   2700,2670,2700 
2700  NCT=2 

IF  <F(6)  -TF)  2710,2720,2710 
2710  PRINT  1330 

STOP  20 
2720  PRINT  1340. TC,TF 

IF  (ND-2)  2750,2730,2750 
2730  PRINT  1010 

PRINT  1350 

PRINT  1360,  (J,KP( J) ,J  =  1  ,NE) 
;      COMPUTE  CRITERIA  FOR  RK  II  STEP  SIZE 

DO  2740  J=1,NE 
2740      AU)  =10.0**KP(  J  )«  1.  0E-04/C  TF-TO  ) 

AC31 )=TF-TO 
:      INITIAL  STEP  SIZE  SET  FOR  RK  II 

F(3)=A(31 )»1.0E-05 

XI49)=F(3) 

GC  TO  2770 
2750  PRINT  1370,  F(3),  TO,  FJ4) 

GO  TO  (2770, 2760), NDT 

PRINT  1370,  F(5),  F(4),  F(6) 

PRINT  1020 

READ  PRINTOUT  CARD 

READ  1380,  (ICHECK(I), 1= 1 ,6 ) , INCPR, ( ICHECK ( J ) , J=7, 1 3 ) 

ICHECKI 14)=8HCATA  LIS 

ICHECKI 15)=4HHELC 

ICHECKI  16)=4HNCNE 

ICHECK ( 17 )= ICHECKI 3 )+ICH£CK{ 10) 

ICHECKI  18)  =  ICHECM3)  +  ICHECK(  12) 

I  CHECK ( 19)  =  ICHECKI 10)  + ICHECKI 12) 

ICHECK(20)=ICHECK( 14) ♦ICHECKI 15) 

ICHECK ( 21 )=ICHECK( 14)+ ICHECKI 16) 

ICHECK I  22)  =  ICHECK I  15 ) +  IC HECK ( 16  ) 

IF  (  ICHECK(  17J-ICHECM20  )  )  2780,2800,2780 
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2780  IF  {  ICHECM18)-[CHECM21  )  )  2790,2800,2790 
2790  IF  {  ICHECKC19)-ICHECM22  )  )  2810,2800,2810 
2800    PRINT    1M00 

STOP    22 
2810    IF     (  ICHECK(16)-ICHECK( 12) )    2320,2940,2820 
2820    IF    UCHECM15)-ICHECM  10)  )    2840  ,  2830  ,  2840 
2830    INCPR=IR1 

GC    TO    2940 

2840  IF  (  ICHECK(14)-tCHECK{3)  )  2841,2050,2841 

2841  PRINT  1390 
STOP  21 

IF  (INCPR)  2880,2841,2880 
IR1=INCPR 

READ  X  VALUES  TC  BE  PRINTED 
READ  1410,  (  IP(K),K=1,8) 

DC    2900    1=1,8 

IF    UPIIII    2900,2890,2900 
2890  NP=I-1 

GC    TC    2910 
CONTINUE 
READ    PRINTOUT    TITLE    CARDS 
READ    1420,    NF,     ( JT ITLEI J ) ,NH , J= 1 ,8 ) 
READ    1420,    NF,     ( JT 1TLE1K ) ,NH ,K=9, 1 6 ) 
;  READ    GRAPH    CARD 

2940  READ    1430,     ICHECM  1 ),  ICHECM  2  ),  NG,  ICHECM  3  ),  ICHECM4  ),  ICHECM  5  ),  NF 
1, (ICHECM I ), 1=6,13) 

ICHECM  14)=8F    SINGLE 

ICHECM  15)=8F    MLLTIPL 

ICHECM 16)=4FHELD 

ICHECM 17)=4FN0NE 

I CHECK ( 18)  =  ICHECM3)+ICHECM6) 

I CHECK C 19)=ICHECM3)  +  ICHECM 10) 

ICHECM20)  =  ICHECK(3)MCHECM 12) 

ICHECM  21  )  =  ICHECM  6)* ICHECM  10) 

ICHECM  22>  =  ICHECM6)+ICHECM  12) 

ICHECM  23)  =  ICHECM 10 ) *  ICHECK ( 12 ) 

ICHECKJ  24  )  =  ICHECM  14)  + ICHECM  15) 

ICHECM  25)  =  ICHECM  14)*  ICHECM  16) 

I  CHECK ( 26 )  =  I CHECK (  14  )  +  ICHECK ( 17 ) 

ICHECM  27  )  =  I  CHECK  (  15  )♦  ICHECK  {  1  6  ) 

ICHECM  28)  =  ICHECM  15)-UCHECM  17) 

ICHECK (29)= ICHECK ( 16)+ ICHECK { 17) 

IF    (  ICHECM  18)-ICHECM24  )  )     2941,2950,2941 

2941  IF     (ICHECK(19)-ICHECK(25) )     2942,2950,2942 

2942  IF    (  ICHECM2C)-ICHECM26)  )    2943,2950,2943 

2943  IF    (  ICHECM21)-ICHECK(27)  )     2944,2950,2944 

2944  IF    (  ICHECM22)-ICHECK(28)  )    2945,2950,2945 

2945  IF    ( ICHECK(23)-1CHECK(29 ) )    2960,2950,2960 

2950  PRINT    1440 
STOP    24 

2951  PRINT    1440 
STOP    25 

2960  IF    UCHECM3)-ICHECK{14)  )    2962,2961,2962 

2961  MCD=1 

GO    TO    3000 

2962  IF    (  ICHECM6)-ICHECK(15)  )    2970,2963,2970 

2963  NG=NH 
MCD  =  2 

GC    TO    3000 
2970    IF    {  ICHECM17)-ICHECM  12  )  >    2980,3130,2980 
2980    IF    (  ICHECM  16)-ICHECM10))    2951,2990,2951 
2990    NG=IR2 

GO    TO    3130 

3000  IR2=NG 

;  READ    PLOT     INCREMENT    CARD 

READ    1460,     ICHECM  1  )  ,  ICHECM2),  INCGR,  ICHECM  3) 

ICHECM4)=8HPL0T    EVE 

IF    (  ICHECM1  J-ICHECM4))     3001,30  10,3001 

3001  PRINT    1510 
STOP    26 
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3010  IHNG)  3020,3011,3020 

3011  PRINT  1500 
STOP  28 

3020  IFCNG-5)  3040, 3040, 3030 
3030  PRINT  1480 

STOP  26 
C      READ  X  VALUES  TO  BE  PLOTTED 
3040  NG2=NG*2 

READ  1470.  (KHtlG(K)  ,NH,  IGCK+1)  , NH, K=l , NG2,2  ) 

IF  (INCGR)  3050,3050,3130 
3050  PRINT  1490 

STOP  27 
3130  PRINT  1520 

IF  (INCPR)  3150,3140,3150 
3140  PRINT  1530 

GO  TO  3180 

3150  IF  (INCPR-1)  3151,3151.3152 

3151  IF  (NP-11  3153,3153,3154 

3152  IF  (NP-1)  3155,3155,3160 

3153  PRINT  1543 
GO  TO  3170 

3154  PRINT  1541 
GO  TO  3170 

3155  PRINT  1542,  INCPR 
GC  TO  3170 

3160  PRINT  1540,  INCPR 

3170  PRINT  1550,  ( IP  (  J  )  ,  JTITL  E  (  J)  ,  JT  ITLE  U  +  08  )  ,  J=  1 ,  NP  ) 

3180  PRINT  1020 

PRINT  1570 

IF  (NGJ  3200,3190,3200 
3190  PRINT  1580 

GC  TO  3480 
3200  GC  TO  (3210,3300  , MOD 
C      SINGLE  GRAPH  SECTION  PRINTOUT 

3210  IF  (NG-1J  3211,3211,3212 

3211  PRINT  1591 
GC  TO  3213 

3212  PRINT  1590. NC 

3213  IF  (INCGR-1)  3214,32  14,3215 

3214  PRINT  1611 
GC  TO  3216 

3215  PRINT  1610, INCGR 

3216  PRINT  1620, IG( 1 ),IG(2) 
IF  (NG~2>  3260,3220,3220 

3220  PRINT  1630,IG<3).IG{4) 

IF  (NG-3)  3260,3230,3230 
3230  PRINT  1640. IC(5),IG(6) 

IF  (NG-4J  3260,3240,3240 
3240  PRINT  1650* IG(7),IG(8) 

IF  <NG-5>  3260,3250,3250 
3250  PRINT  1660, IG(9),IG( 10) 
3260  DC  3270  J=1,NG2 

IF  UG(J)-5C)  3270,3290,3270 
3270      CONTINUE 

GO  TO  3480 
3290  PRINT  1560 

GO  TO  3480 

MULTIPLE  GRAPH  SECTICN  PRINTOUT 

IF  (NG-1)  3301,3301,3302 

PRINT  1601 

GC  TO  3303 

3302  PRINT  1600, NG 

3303  IF  (INCGR-1)  3304,3304,3305 

3304  PRINT  1611 
GO  TO  3306 

3305  PRINT  1610, INCGR 

3306  PRINT  1670,  IGC1).IG{2) 
IF  (NG-2J  3260,3310,3310 

3310  PRINT  1680,  IG(3),IG(4) 
IF  (NG-3)  3260,3320,3320 
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3320 
3330 
3340 
3480 

3500 
3510 

3520 
3521 
3522 

3523 


3524 


3525 

3530 
35U0 
3550 
3560 
3570 
3580 


3590 
3600 
3601 

3610 


3620 

3630 

3640 

3650 
3660 
3670 
3671 

3680 
3730 


PRIN 
IF  ( 
PRIN 
IF  ( 
PRIN 
GC  T 
PRIN 
PRIN 
PRIN 
IF  ( 
PRIN 
PRIN 
T  =  T0 
DT  =  F 
DO  3 

GO  T 
DO  3 
D(K) 

MYCK 

C(50 

RETU 

SET 

MY{  1 

MA=K 

MB  =  y 

MOM 

MD=M 

ME  =  M 

MF  =  M 

MG=M 

MH=M 

MI  =  K 

LINE 

NCPT 

NUMP 

ENTR 

CC5C 

RETU 

IF  ( 

IF 

IF 

IF 

IF 

PRIN 

NPAG 

PRIN 

PRIN 

PRIN 

GO  T 

C(50 

RETU 

LINE 

DO  3 


T  1690,  IG(5),IG(6) 
NG-4)  3260, 3330, 3330 
T  1700,  IG(7),IG(8) 
NG-5)  3260,33i40,33UO 
1710,  IG!9),IG110) 


3260 
1C00 
1730, 
1740, 
I NC PR) 
T  1750, 
T  1010 


UTITLE(J),  J*l,  5) 
NRC 
3500,3510,3500 
(JTITLEt J),J=1,16) 


35 
35 


PRIN 
IF  ( 
IF  i 
GG  T 
C(50 
RETU 
GRAP 
DC  3 


S3) 

520  J=1,N 

X(J)=XC(J 

0  43521,3 
522  K=l,3 
=0. 

1  =  1 

)  =  5. 

RN 

INDEXES  F 

)  =  2 

YI2) 

Y(3) 

Y(U) 

YC5) 

Y(6) 

Y(7) 

Yf8) 

Y(9) 

Y{  10) 

S  =  0 

S  =  0 

TS=0 

Y  POINT  F 

)=6. 

RN 

INCPR) 

NGPTS) 

XKODFCNOP 

XMODF(NOP 

XNODFCNOP 

T  1000 

E=NPAGE+1 

T  1741,  N 

T  1750»  i 

T  1010 

C  {3601 

)=2.0 

RN 

S=LINES+1 

640  J=1,N 

IF  (  IPU  ) 

PR(J)=T 

GO  TO  36M 

IPJ=IPf J 

PR«  J)=X( 

CONTINUE 

T  1760e 

NG1  3660 

XMODF«NGP 

0  C  3671, 2 

)=3.0 

RN 

H  POINTS 

750  J=1,K 

IF  <IG(J) 

GRC J)=T 


E 
) 
524), KZ 

0 


CR  SPECIAL  FUNCTIONS 


ROM  RK  I  OR  RK  II 


,3 


PACE, NRC 

JTITLE( J),J=1, 16) 

61C),NA 


-5C)    3630,3620,3630 

0 
i 
IPJ) 


PR( J),J=1,NP) 

3850,3660 

TS.INCGR))  3850,3670,3850 

68C),MA 


STORAGE 

G2 

-50)    37M0, 3730,3740 
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GC    TC    375C 
37i<0  IGJMG(J) 

GRi J)=X{  ICJ) 
3750    CONTINUE 
3840    NUMPTS  =  NUN,PTS+1 

X1(NUMPTS)=GR( 1 ) 

Y1(NUMPTS)=GR(21 

X2CNUMPTS)=GRC3) 

Y2CNUMPTS)=GP(4) 

X3(NUMPTS)-GRC5J 

Y3fNUMPT$)=GPC6) 

X4CNUMPTS )=GRC7) 

Y4(NUMPTS)=GR(8) 

X5tNUMPTS)=GR(9) 

Y5INUMPTS)=GRC 1C) 

3850  N0PTS=NCPTS+1 

IF  CLINES-30C)  3852,3851,3851 

3851  PRINT  1770 
GC  TC  8000 

3852  IF  (T-l.E+04)  3870,3860,3860 
3860  PRINT  1780 

GC  TC  8000 
3870  IF  (NCPTS-IOCOC)  3890, 388C, 3880 
3880  PRINT  1790 

GC  TC  8000 
3890  IF  (NUMPTS-9C0)  3910,3900,3900 
3900  PRINT  1800 

GC  TO  8000 

IF  (T-IF)  3930,3920,3920 

PRINT  1810 

EXIT  TO  GRAPH  PLOTTING  ROUTINE 

GC  TO  8000 
3930  DC  3950  J=1,49 

IF  ?ABSF(X( J))-l.E+04)  3950  ,  3940  ,  3940 
3940      PRINT  182CJ 

GO  TO  80C0 
3950      CONTINUE 
C      SPECIAL  FUNCTION  SECTION 

GC  TO  14160, 3960) ,MA 
3960  GO  TC  ( 4090,3970,  MR 
3970  GO  TO  f  4010,3980, fC 
C      DELAY1 
3980  IF  (MYM1)-5C0)  4000,3990,3990 
3990  MYC 1 1)=0 
4000  MYC 1 1 )=MY  C 1 1  )  +  l 

MX=MYC 1 1) 

VCMX )=D( 3) 

W«MX)=DC  1  J 

MYC  14)=MYC 15) 

MYC 12)=MYt 13) 
4010    GO    TO    { 4050,14020,  MO 
C  DELAY2 

4020    IF     (MYC 16I-1CO0)    4040,40  30,4030 
4030    MYC 16)=500 
4040    MY* 16)=MY^ 16  1+1 

MX=MY« 16) 

VIHX)*0(7J 

W!MX)=Df 5) 

MYC 19)^MYC20) 

MY! 17)=MYC 18) 
4050  GC  TC  «4090,4C6C),ME 
C      DELAY3 
4060  IF  (MY<21 H-1500)  M08C, 4070, 4070 
4070  MYC 2 11=1000 
4080  MYC21)=MYI21 H 1 

MX=MYC21 ) 

VCMX) =0(1 1  I 

W(MX)=D(9) 

MYC24)=MYC 25) 

MYC22)=MYC23) 
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U090 

14  1  00 

4110 
'4120 

4130 
i+IUO 


4150 
4160 


4170 


GO 

REL 

t)l  1 

DC  1 

D{  1 

GC 

REL 

DM 

DC2 

DC2 

GC 

RKL 

DC2 

MJ  = 

GC 

D«2 

TAD 

DC 


EXI 
GC 


TC  UllC,inCC),MF 

30W  1,2  am;  3 

4j=0' 13) 

6>=D( 15) 

B)  =  D( 175 

TC  U130,M12C),MG 

2H  1,2  ANC  3 

9)=D(20) 

1  )=D(22) 

3)=D(24) 

TO  ft4l50,lil4C},MH 

ASH 

5)=D(26) 

MYC26) 

TC  (  4160, 4150), MJ 

71=0(26) 

=  T 

4170  L=1,NE 

XA{L)=X(L) 

XBIL)=Xf L) 
T  TO  RK  I  CR  RK  I  I 
TC  I  5000,7000  ,ND 


BEGINING  CF  RUNGA  KUTTA  I 


5000  DC  5030  1=1,4 

T=TAD+CT? I )»0T 
DC  5010  J=1,NE 

XU)=XA?  J)+CT(  I)*AK(I-1,J) 
C(50)=4.C 
RETURN 
DC  5030  J=1,NE 

AKI  I,  J)*DT»XDGTC J) 
DC  5040  vJ=1  ,NE 

XA(J}=XA(J)"MAKf1,J)+2.*AK(2,J)*2.*AK(3,J)*AK«4,J))*0.  1666666667 
GO  TO  C6000, 7010,7040, 7060) ,IND1 


5010 


5020 
5030 

5040 


;      END  CF  RUNG  KCTrA  I 

6000  IF  (NCT-1)  6C20, 6010, 6020 
6010  DT=F(3) 

GC  TC  6050 
6020  IF  (NDT-2)  6C40 , 6030 , 604 0 
6030  IF  IT-FU1)  6010,6040,6040 
6040  DT=F(5) 

X(49)=0T 
6050  T=TAD+DT 

DC  6060  J=1,NE 
6060      XU)=XAU) 
;      LEAVE  RK  I  TC  PRINT  AND/CR 

GC  TC  13530,3525), MA 


PLCT  POINT 


BEGINNING  CF  RUNGA  KUTTA  II 

THIS  SUBROUTINE  IS  BASED  ON  THE  METHOD  DESCRIBED  IN  CCMP.  CF 
THE  A.C.M.,  JUNE  1960  (PP.  355  -360).  IT  CHECKS  THE 
TRUNCATION  ERROR  AT  EACH  STEP  CF  THE  INTEGRATION  AND  ADJUSTS 
THE  STEP  SIZE  ACCORDINGLY,  THE  ACTUAL  RUNGE-KUTTA  INTEGRATION 
IS  PERFC90t4  BY  RKUTTA,  WHICH  IS  CONTAINED  IN  THIS 
PROGRAM  (5CC0  SERIES).   THE  ARGUMENTS  ARE, 

NE  =  NUMBER  CF  (FIRST  ORDER)  EQUATIONS.   MAXIMUM  NE=30 
T  =  TIME  AT  START  OF  INTEGRATION  STEP  (UPDATED  BY  THIS 
SUBROUTINE  AFTER  THE  COMPLETION  OF  EACH  STEP). 
DT  =  STEP  SIZE  AT  START  CF  EACH  STEP  (ALSC  UPDATED). 
XU)  -  THE  N  DEPENDENT  VARIABLES  (ALSO  UPDATED). 
A(I)  =  THE  SPECIFIED  ALLOWABLE  ERROR  PER  UNIT  TIME  FOR  EACH 
CF  THE  DEPENDENT  VARIABLES. 
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C  A?31HS  USED  TO  ENTER  I  HE  TOTAL  TIME. 

C  NOTE  THAT  IF  NECESSARY  DT  IS  REDUCED  FROM  ThE  VALUE  STATED 

C  IN  THE  ARGUMENT  UNTIL  THE  SPECIFIED  ACCURACY  HAS  BEEN 

C  ACFIEVEC. 

7000  PREVDT  =  f:r 

TS=T 

H  =  CT 

DT=2.G*h 

IND 1=2 

GC  TO  5000 
7010  LC  7020  1=1, NE 

XRi  I)=XA{ ! ) 
7020  XAU)=XBU) 
7030  DT=H 

I  NCI =3 

GC  TC  5000 
7040  DC  7050  1=1, NE 
7050      XRRU)=XMI) 

TAD=TS*DT 

INDI=4 

GC    TC    5000 
7060    DO    7070     1=1, NE 

XS(I)  =  XAU) 
7070  Xf I )=XBI  I) 

U2=0.03 

DC    7100     1=1, NE 

E21=SXR$  D-aSIIS 1*0.0666  6666  667 

E2lR  =  niMF(AHSF(E21  )t  ABSFUARSF(XS<  I  >  MABSFi  XRU  ))-1.99« 
1  ABSFiXC  I)))*1.0E-09)  ) 

C  THIS    TENDS    TC    PREVENT    ROUND    OFF    ERROR    FROM    TAKING    CONTROL 

U=    E21R    /    (A(|)«6.0*HJ 
IF     SU-U2)    7100,7100, 7C90 
7090  U2=U 

7100  XSf  H  =  XSf  D-E21 

IF    IU2-1cCJ    7150,7110,71 1C 
7110    IF    (H-A131 >*1.0E-09J    7120,7120,7130 
7120    DT=AJ31 )*1.0E-09 
C  THIS    SETS    THE    MINIMUM    STEP    SIZE    TO     l.OE-09    TIMES    THE    TCTAL    TIKE 

GC    TC    7200 
7130    DC    7  HO    I  =  1,NE 

XAC I)=XBC I  I 
7140  XRm=XRfHI) 

H=0.5*H 

TAD=TS 

GC  TO  7030 

THIS  RECYCLES  THE  INTEGRATION  IF  THE  TRUNCATION  ERRCR  IS  EXCESSIVE 

IF  IU2-0.031)  7160,7170,7170 

DT=2oC*H 

GC  TC  7180 
7170  CT*SCRTFCSCRTF(C.5/U2) )*H 
7180  IF  10T-AI31)*C.C01)  7200,7200,7190 
7190  DT=A«31 }*OoQCl 

GC  TC  7230 

THIS  SETS  THE  MAXIMUM  STEP  SIZE  TO  1.0E-03  TIMES  THE  TCTAL  TIME 

IF  (DT-PREVDT)  7230,7210,7230 

ICCUNT=ICCUNT+1 

IF  UCOUNT-3)  7230,7220,7220 
7220  ICOUNT=0 

DT=AC31 )«l,0E-05 
;      THE  ABOVE  IS  INTENDED  TO  CLEAR  A  STEP  SIZE  JAM 
7230  T=TS+2.0*H 

XCU9)=DT 

DC  72U0  1=1, NE 
XI  I )  =  XS?  I) 

THIS  UPDATES  T  AND  XII). 

DT  IS  UPDATED  BY  STATEMENTS  7120,  7160,  7170  AND  7190 

END  CF  RUNG A  KUTTA  II 

LEAVE  RK  II  TC  PRINT  AND/CR  PLOT  POINT 

GC  TC  {3530,3525), MA 
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C      GRAPH  PLOTTING  ROUTINE  I  6CCC  SERIES) 

8000  IF  SNC)  8010,8320.8010 

8010  GC  TC  S 8020, 8C3C,80U0, 8050,8060, 8070,8090, 81CC,81 10) ,NRC 

8020  ITITLEC6)=BH   RUN  1 

GO  TC  8120 
8030  ITITLEI6)=8H   RUN  2 

GC  TC  8120 
8040  ITITLEC6)=8H   RUN  3 

GO  TC  8120 
8050  ITITLE«6)=8H   RUN  k 

GO  TC  8120 
8060  ITJTLEC6)=3H   RUN  5 

GO  TC  8120 
8070  ITITLE?6)=8H   RUN  6 

GC  TO  8120 
8090  ITITLE?6?=8H   RUN  7 

GC  TO  8120 
8100  ITITLEI68=8H   RUN  8 

GC  TO  8120 
8110  IT!TLEC6)=8H   RUN  9 

8120  SFX=0. 
SFY=0. 
MINCFFX=0 
MINCFFY=0 
LABELN0=11 
MCDE=0 

N2  =  0 
N1=0 

GO    TO    !8121f  819O.M0C 
C  SINGLE    GRAPH    SECTION 

8121  MCDCURV=0 
LABEL=4H 
KCL=1 

8130  ITITLE « 7J=8H  GRAPH  A 

8131  CALL  GRAPH2  \ NUVPTS, X 1 , Y 1 ,8 , MODCURV , LABEL , ITI TLE , SFX ,SFY, MlNOFFX, 
1  MI\0FFY,LABELN0,M0CE,N1 ,N2) 

GO  TO  (8132. 6320, 8220),KCL 

8132  IF  (NG-2)  8 1 80,8 140, 81 40 

8140  ITITLE(7)=8H  CRAPH  B 

8141  CALL  GRAPH2  CNU^PTS, X2 ,Y2 , 8, MODCURV , LABEL, IT ITL E, SFX , SFY,  MlNOFFX , 
1  MINCFFY,LABELNC,M0CE,N1,N2) 

GO  TO  C8142, 8320, 8250), KCL 

8142  IF  (NG-3)  8180,8150,8150 

8150  ITITLE(7)=8H  GRAPH  C 

8151  CALL  CRAPH2  ( NUMPTS, X3 ,Y3, 8, MCDCURV , L ABC L, IT  I TLE , SFX ,SF Y,  MI NCFFX , 
1  M  INOFFY  ,  LAB  ELNC  ,MOi)E ,  N  1  ,  N2  ) 

GC  TO  18152,8320,8280) , KCL 

8152  IF    CNG-4)     8180,8160,8160 

8160  ITITLEC7?=8H  GRAPH  0 

8161  CALL  GRAPH2  ( MJMPTS , X4, YU ,8 , MODCURV, LABEL , IT  I TLE, SFX, SFY, M INOFFX, 
1  MI?\CFFY,LABELN0,M0CE,N1,\'2) 

GO  TO  «8162,832C«8310) ,KOL 

8162  IF  (NG-5)  8180,0  170,8170 

8170  ITITLEC7i=8H  GRAPH  E 

8171  CALL  GRAPH2  iNUMPTS , X5,Y5 ,8 , MODCURV. LABEL, IT  I TLE, SFX , SFY, M INOFF X, 
1  MIF\0FFY,LABELNC,M0CE,N1  ,N2) 

8180  GC  TO  8320 
C      MULTIPLE  GRAPH  SECTION 
8190  ITITLE(7)=8H  GRAPH  M 

LABEL=4H  A 

IF  CNG-1)  82CC, 8200, 8210 
8200  MCDCURV=0 
K0L=2 

GO  TC  8131 
8210  MCDCURV=1 

KCL  =  3 
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GO  TG  8131 
8220  MCDCURV=2 

LABEt=UH  R 

IF  (NG-2)  8230, 8230, 82U0 
8230  M0DCURV=3 

KCL=2 

GC  TO  81U1 
82^0  K0L=3 

GC  TG  81U1 
8250  LABEL=UH  C 

IF  CNG-3)  8260,8260,8270 
8260  MCDCURV=3 

KCL  =  2 

GC  TC  8151 
8270  K0L=3 

GO  TO  8151 
8280  LABEL=4H  D 

IF  (NG-U)  8250,8290,8300 
8290  MCDCURV=3 

KCL  =  2 

GC  TC  8161 
8300  KCL=3 

GO  TC  8161 
8310  MCDCURV=3 

LABEl=4H  E 

GC  TC  8171 
8320  PRINT  1000 

IF  INRC-NRJ  215C, 8330, 8330 
8330  PRINT  1830 

STOP 

END 

END 
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AVAILABLE  FUNCTION  SUBPROGRAMS 


FUNCTION  DELAYl  (DELT,VAR,T) 

DIMENSION  V(  1S00),W{ I  500 } , MY ( 30 ) , D{ 4  0 ) 

COMMON  V,W».MY,D 

MA  =  MYC 1  ) 

GO  TO  ( 1,2), MA 

1  D(!)*VAR 
DC2>=VAR 
D(3)=T 
MYC2)=2 
MYI3)=2 
MY?4)=2 
MY( 1  1)=0 
RETURN 

2  DC 1 )=VAR 
D(3)=T 
MC=MY{ 12) 

GO  10  (3,9),MC 

3  MYC 13)=1 

IF  (CELT-T)  8tU,U 
«♦  IF  (MYM1)-5C0)  7,5,5 

5  PRINT  6 

6  FORMAT  C43H  PERIOD  OF  DELAYl  EXCEEDS  AVAILABLE  MEMORY.) 
STOP  31 

7  DELAY1=D(2) 
RETURN 

8  MYi 13)=2 

9  ME  =  MYHU) 
U(4)=T-DELT 
IC0UNT=1 

10  DO  12  K=ME,5C0 

IF  IV!K)-C(U))  12,12,11 

11  MF=K 

MYf 15)=K 
GO  TO  15 

12  CONTINUE 
ME=1 
ICOUNT=ICnUNT+l 

IF  UCOUNT-3)  1C,  13,  13 

13  PRINT  1U 

1U  FORMAT  C63H  CELAY1  MALFUNCTIONS  CHECK  THAT  DELT  EXCEEDS  INITIAL  ST 
1EP  SIZE.) 
STOP  32 

15  IF  (MF-1)  16,16,17 

16  11=500 
12=1 

GO  TO  18 

17  Il=MF~l 
I2=MF 

18  DELAY  1=M  II  )  *{W(  I2)-W{  II  )  )  *  (  DU  )~V(  I  1)  I  /  i  V  C  I  2  )-VU  1  )  ) 
END 

END 


FUNCTION  DELAY2  { DEL T, VAR , T ) 

DIMENSION  V(  150C),W( 1 500  )  ,MY f 30 ) fDC UO ) 

COMMON  V,W,MY,D 

MA  =  MYC 1  ) 

GO  TC  i  1,2), MA 

D(5)=VAR 

D(6)=VAR 
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DC7)=T 

MYf  2)  =2 

MY* 3) -2 

MYJ5)=2 

MY( 16)=500 

MY! 1 9 ) =50 1 

RETURN 

D(5)=VAR 

I)  (  7  )  =  T 

MC=MY( 17) 

CO  TC  C3,9),MC 

MYS 18)  =  1 

IF  (OELT-T)  8,4,4 

IF  (MY*  16)-1C00)  7,5,5 

PRINT  6 

FORMAT  S43H  PERIOD  OF  DELAY2  EXCEEDS  AVAILA3LE  MEMORY.) 

STOP  33 

DELAY2=D(6) 

RETURN 

MYf 18)-2 

ME=MY? 19) 

DC8)=T-DELT 

ICCUNT=1 

DO  12  K=ME,1CC0 

IF  (V(K)-ni8)>  12,12,11 

MF  =  K 

MYC20)=K 

GO  TO  15 
12      CONTINUE 
ME=501 

IC0UNT*IC0UNT+1 
IF  UCCUNT-3)  10,13,13 
PRINT  14 

FORMAT  J63H  CELAY2  MALFUNCTION.  CHECK  THAT  DELT  EXCEEDS  INITIAL  ST 
1EP  SIZE.) 
STOP  34 

IF  (MF-  501)  16,16,17 
11=1000 
12=501 
GO  TO  18 
I1=MF-1 
I2  =  MF 

DELAY2  =  V««  II  )*(W<  I2)-M  II  )  )*(D(8)-V(  II  )  )/(V(I2J-V(I1  )  > 
END 
END 


10 
11 


13 
14 


15 
16 


17 
18 


FUNCTION  0ELAY3  (DELT,VAR,T) 

DIMENSION  V«  1500) VW( 1500 )  ,MY(30) ,D 

COMMON  V,W,M.Y,D 

MA=MYU  ) 

GO  TO  U,2),MA 

D(9)=VAR 

D(10)=VAR 

D(11)=T 

MY(2)=2 

MY(3)=2 

MYI6)=2 

MY«21 )=1000 

MYC24S=1001 

RETURN 

DI9)=VAR 

DC  1  1  )  =  T 

MC=MY«22) 

GO  TO  £3,9) , 

MY*23)=1 

IF  (DELT-T) 

IF  (MY(21)-1 

PRINT  6 


40) 


MC 

8,4,4 

500)  7,5,5 
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6  FORMAT  U3H  P[RI(0  OF  DELAYS  EXCEEDS  AVAILABLE  MEMORY*  I 
STOP  35 

7  DELAY3~0( TO) 
RETURN 

8  MY* 23)  =  2 

9  MC=vy«2U) 

Df 12)=T-DELT 
[COUNT*! 

10  DC  12  K  =  ME, 15C0 

IF  (V{K)-0( 125  J  12,12,1  1 

11  MF  =  K 
MYS25)=K 
GO  TC  15 

12  CONTINUE 
ME=10G1 
ICCUNT=ICOUNT*  1 

IF  {  ICOUNT-3)  10,13,  13 

13  PRINT  14 

1U  FORMAT  (63H  LELAY3  MALFUNCTION.  CHECK  THAT  DELT  EXCEEDS  INITIAL  ST 
1EP  SIZE.) 
STOP  36 

15  IF  (MF-1001)  16,16,17 

16  11=1500 
12=1001 
GO  TO  18 

17  I1=MF-1 
I2  =  MF 

18  DELAY3=WC  1 1  I  ♦  C W (  12  )-W(  II  )  ) « ( 0( 1 2 )-V ( 1 1 ) ) / { V(  12 ) -V (  II ) ) 
END 

END 


FUNCTION  OZONE  ( X  INPUT, H ALFO, GA I N ) 

IF  tABSFCX INPUT l-HALFD)  1,2,2 

DZONE=0. 

RETURN 

0UT=GA1N*IABSFC  < INPUT )-H ALFO ) 

DZCNE=SIGNF(CLT,X[NPUT) 

END 

END 


FUNCTION  REL2  { X  INPUT, CU TPUT ) 
REL2=SIGNFC0LTPUT, XINPUT ) 
END 
END 


FUNCTION  REL2H1  ( X  INPUT, HALFD, OUTPUT) 

DIMENSION  MY(30),fMU0) 

COMMON  MY,D 

MA=MYM  ) 

GO  TC  « 1,2), MA 

1  MY«2)=2 
MY«8)=2 
DU9)=0UTPUT 

2  IF  ?  DC  1 9  S )  5,3,8 

3  PRINT  4 

4  FORMAT  (20H  PEL2H1  MALFUNCTION.) 
STOP  37 

5  IF  (XINPUT-HALFC)  6,7,7 

6  REL2hl=D( 19) 
DC20)=D(19) 
RETURN 

7  REL2H1=0UTPUT 
D(20)=0UTPUT 
RETURN 
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3 
9 


10 


IF  (X  INPUT  •H/'LFL  ) 

REL2Hl=Ct 19) 

0(20=0(191 

RETURN 

REL2H1=-0UTPUT 

Df2C)=-CUTPUl 

END 

END 


10.  10,9 
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FUNCTION  REL2H2 

DIMENSION  MY(30) 

COMMON  MY,D 

MA=MY! 1  J 

GO  TC  !1,2),NA 

MYC2)=2 

MY(8)=2 

DC  21  )=OUTPUT 

IF  10(21))  5,3,0 

PRINT  4 

FORMAT  I20H  REL2H2 

STOP  38 

IF  (XINPUT-HALFU)  6,7,7 

REL2H2=Di21) 

0! 22) =0  121 3 

RETURN 

REL2H2=0UFPUT 

0(22)=0UTPUT 

RETURN 

IF  (  X  INPUT  *H At. FO)  10,10,9 

REL2H2=U(21) 

D(22)=DC21 5 

RETURN 

REL2H2=-OUTPUT 

D(22)=-OUTPUT 

END 

END 


IX1NPUT, HALFD, OUTPUT) 
,  D  (  4  0  ) 


MALFUNCTION, I 
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FUNCTION  REL2H3 

DIMENSION  MY(30) 

COMMON  MY,D 

MA=MY£ 1  ) 

GO  TO  f 1,2),^A 

MYf2)=2 

MYC8)=2 

DI23)=0UTPUT 

IF  (0(231)  5,3,8 

PRINT  1* 

FORMAT  (20H  PEL2H3 

STOP  39 

IF  (XINPUT-HALFU)  6,7,7 

REL2H3=D(23) 

D(24)=D(23) 

RETURN 

REL2H3=CUTPUT 

0(24)=CUTPUT 

RETURN 

IF  (XINPUT+HALFC)  10,10,9 

REL2H3=UI23) 

D(24)*D(23) 

RETURN 

REL2H3=-OUTPLT 

D(24)=-GUTPUT 

END 

END 


(XINPUT,HALFD, OUTPUT) 
,0(40) 


MALFUNCTION. ) 
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FUNCTION  REL3  I XINPUT, OUTPUT  J 
IF  {X  INPUT?  2,1,2 

1  REL3  =  CO 
RET  U R N 

2  REL3=SIGNF JCUTPUTtXINPUT ) 
END 

END 


FUNCTION    REL3D     (  X  INPUT, HALFli, OUTPUT) 
IF    i ABSFC XINPUT  )-HALFO)     1,2,2 

1  REL3C  =0.0 
RETURN 

2  RE13D=SIGNF{CUTPUT*X  INPUT ) 
END 

END 

FUNCTION  REL3DH1  ( XINPUT, DROPOUT, PULL  IN, OUTPUT ) 

DIMENSION  MY  130)  ,DC40) 

COMMON  MY,D 

MA=MYi  1  J 

GC  TC  H,2),MA 

1  MY£2)=2 
MY£7)=2 

2  IF  lAPSFUINPUD-DRCPOUT)  3,4,4 

3  REL3DM=0. 
DC  13  1=0. 
RETURN 

4  IF  (APSFU  INPUT  J-PULl  IN)  6,5,5 

5  Df  1 3  )  =  SIGNFlCLTPUT, XINPUT) 
REL3l)H=U( 13) 

RETURN 

6  REL3CHl=nnU) 

DC  i3)=nnun 

END 
END 


FUNCTION  REL3DH2  ! XINPUT .DROPOUT, PULLI Nt OUTPUT ) 

DIMENSION  MY?30),D£40) 

COMMON  MY,D 

MA=MY? 1  ) 

GC  TO  H,2)  ,MA 

1  MYf2>=2 
MY(7)=2 

2  IF  I APSFf .  XINPUT)-OROPOUT )  3,4,4 

3  REL3DH2=0. 
D!15)=0. 
RETURN 

4  IF  (APSFCXINPUT)-PULLIN)  6,5,5 

5  D  U  5  )  =S I GNF i OUTPUT, X INPU T ) 
REL3DH2  =  IM 15) 

RETURN 

6  REL3DF;2  =  DZ  16) 
D(15)=0{16) 
END 

END 


FUNCTION  REL3DH3  I XI NPUT , DROPOUT , PULL  IN, OUTPUT ) 

DIMENSION  MY {30) ,0540) 

COMMON  MY,D 

MA=MYK 1 ) 

GC  TO  {  1,2), MA 

*Yf 25=2 
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IF 

RE 
DC 
RE 
IE 
D( 
RE 
RE 
RE 
Df 
EN 
EN 


(7)  = 

(AB 

L3DH 

17)  = 

TURN 
CAP 
17)  = 
L3CH 
TURN 
L3DH 
17)  = 
D 
D 


SFCXINPUTJ-ORCPCUT )  3,U,4 

3  =  0. 

0. 

SF(XINFUT)-PULLIN)     6,5,5 
SIGNF(CUTPUT,XINPUT) 
3=0 f 17) 

3=Oi 18) 
DM8) 


FUNCTION  SATAN*P  { X  INPUT, CUTCFF, GAIN ) 

IF  (GAIN«ABSF(XINPUT)-CUTCFF)  1,2,2 

SATAMP*GAIN*XINPUT 

RETURN 

SATAMP*SIGNF(CUTOFF, XINPUT) 

END 

END 
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